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1  Introduction 


Since  the  first  space  mission  in  1957  (Sputnik  1),  artificial  objects  of  different  size  appeared 
in  orbits  around  the  Earth  .  Nowadays,  the  number  of  such  objects  is  estimated  to  be  over 
2  millions.  More  than  300000  of  them  have  a  size  greater  than  10  cm  although  only  5%  of 
them  are  catalogued.  The  observation  and  tracking  of  Resident  Space  Objects  (RSOs)  has 
become  a  crucial  task  in  launch  planning  of  new  satellites,  collisions-avoidance  operations  and 
in  general  to  ensure  the  safety  of  operational  satellites.  As  shown  in  Fig.  1  the  number  of  space 
debris  has  been  constantly  increasing.  Particularly,  in  2007  and  in  2009,  two  events  (Fengyun 
anti- satellite  test  and  Iridium-Cosmos  collision)  dramatically  increased  the  number  of  debris 
in  Low  Earth  Orbit  (LEO). 
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Figure  1:  Distribution  of  space  debris 


Fig.  2  (figure  from  [1])  also  shows  that  RSOs  are  more  concentrated  in  low  orbits  (altitude 
below  2000  km)  and  for  inclination  angles  near  90  degrees.  As  a  consequence  of  these  con¬ 
siderations,  tracking  and  prediction  of  RSO’s  orbits  represent  key  tasks  in  Space  Surveillance 
and  Tracking  (SST)  operations. 
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■  Satellites 

■  Rocket  Bodies 


Figure  2:  Distribution  ofRSO  per  orbit  inclination  and  per  altitude  respectively 

One  of  the  issues  related  to  RSOs  tracking  and  classification  regards  data  and  track  asso¬ 
ciation,  which  is  the  problem  of  identifying  and  associating  data  that  is  generated  by  the  same 
object.  A  possible  solution  to  this  problem  consists  of  an  extension  of  the  number  of  parame¬ 
ters  (or  features)  that  are  used  to  characterize  an  object  and  use  them  to  discriminate  it  among 
others. 

These  features  can  be  estimated  using  both  Optical  or  Radar  sensors,  which  typically  per¬ 
form  differently  depending  on  the  operational  conditions  (illumination,  distance,  etc).  Optical 
sensors  produce  a  more  accurate  estimation  of  the  object  angular  position  also  with  objects 
that  are  very  distant  from  the  sensor,  but  are  affected  by  weather  conditions  (fog,  clouds,  rain, 
etc).  Radar  sensors,  instead,  provide  accurate  range  and  range-rate  measurements  and  are 
robust  with  respect  to  weather  conditions,  although  they  are  limited  in  terms  of  object  dis¬ 
tance.  Coherent  radar  have  also  the  ability  to  exploit  the  phase  information  to  estimate  object’s 
Doppler  signatures  and  consequently  estimate  parameters  such  as  their  radial  velocity,  rotation 
speed,  etc.  The  present  technical  report  presents  innovative  algorithms  that  exploits  coher¬ 
ent  radar  signals  to  accurately  estimate  RSO  features  to  uniquely  discriminate  an  object  from 
others.  Effective  radar  signal  processing  based  on  time-frequency  analysis  has  been  proposed 
to  provide  the  basis  for  estimators  to  be  developed  that  ensure  both  accuracy  and  robustness 
against  noise.  Optical  measurement-based  algorithms  are  also  reviewed  that  provide  estimates 
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of  object  parameters  with  the  aim  of  jointly  using  radar  and  optical  features  to  improve  target 
discrimination  performances  with  respect  to  radar-only  or  optical-only  measurements.  Results 
are  presented  that  prove  the  effectiveness  of  the  algorithm  proposed.  This  technical  report  is 
organised  as  follows: 

•  In  Section  2  Time-Frequency  transforms  are  discussed  and  parameters  of  interests  are 
identified. 

•  In  Section  3  the  Radar  signal  model  adopted,  the  signal  generation  process  and  different 
models  of  RSOs  are  introduced. 

•  Section  4  describes  the  algorithms  that  have  been  developed  for  feature  extraction  from 
Radar  data.  Particularly  a  manual  and  an  automatic  algorithm  are  presented  and  dis¬ 
cussed. 

•  In  section  5  elements  of  photometry  are  reviewed  and  a  model  for  optical  features  ex¬ 
traction  is  presented  and  applied  to  the  classes  identified  in  Section  3. 

•  Section  6  concerns  the  proposed  classification  algorithms,  specifically  SVM  and  Bayesian. 

•  In  Section  7,  the  algorithm  performances  are  shown  and  discussed  in  the  case  of  radar- 
only  measurements,  optical-only  measurements  and  after  fusing  information  from  both 
radar  and  optical  measurements. 

•  Section  8  provides  the  conclusions  and  future  developments  regarding  RSOs  orbit  deter¬ 
mination  and  classification. 


DISTRIBUTION  A.  Approved  for  public  release:  distribution  unlimited. 


6 


LINEAR  AND  NON-LINEAR  TIME-FREQUENCY  ANALYSIS 
FOR  PARAMETER  ESTIMATION  OF  SPACE  OBJECTS 


2  Time  Frequency  Transforms  based  RSO  features  extrac¬ 
tion 

2.1  Time  Frequency  Transforms 

The  echo  signals  generated  by  rotating  objects  have  spectral  contents  that  change  with  time 
(non-stationary  signals).  Due  to  lack  of  localized  time  information,  the  widely  used  Fourier 
Transform  (FT)  cannot  provide  time- varying  frequency  modulation  information.  A  Joint  Time- 
Frequency  Analysis  (JTFA)  that  provides  localized  time-dependent  frequency  information  is 
needed  for  extracting  time-varying  motion  dynamic  features.  Time-frequency  Distributions 
(TFDs)  includes  linear  and  bilinear  transforms,  such  as  the  Short  Time  Fourier  Transform 
(STFT)  and  Wigner-Ville  distribution  (WVD)  respectively. 

2.1.1  The  STFT 

The  most  standard  approach  to  analyze  a  signal  with  time-varying  frequency  content  is  to 
split  the  time-domain  signal  into  many  segments,  and  then  take  the  Fourier  transform  of  each 
segment.  This  is  known  as  the  STFT  operation  and  is  defined  as 

STFT(t,u)  —  j  s(t')w(t' —  t)  exp{—jujt'}dt'  (1) 

This  operation  differs  from  the  Fourier  transform  only  by  the  presence  of  a  window  function 
w(t).  As  the  name  implies,  the  STFT  is  generated  by  taking  the  Fourier  transform  of  smaller 
durations  of  the  original  data.  Alternatively,  we  can  interpret  the  STFT  as  the  projection  of 
the  function  s  (t’)  onto  a  set  of  bases  w*(t’  - 1 )  exp  jcut’  with  parameters  t  and  c o  .  Since  the 
bases  are  no  longer  of  infinite  extent  in  time,  it  is  possible  to  monitor  how  the  signal  frequency 
spectrum  varies  as  a  function  of  time.  This  is  accomplished  by  the  translation  of  the  window 
as  a  function  of  time  t  ,  resulting  in  a  2D  joint  time-frequency  representation  STFT  (t,cu)  of 
the  original  time  signal.  The  well-known  spectrogram  defined  as  the  square  modulus  of  the 
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STFT  is  a  well-known  tool  for  time-frequency  analysis.  With  a  time-limited  window  function, 
the  resolution  of  the  STFT  is  determined  by  the  window  size.  There  is  a  trade-off  between  the 
time  resolution  and  the  frequency  resolution.  A  larger  window  has  higher  frequency  resolution 
but  a  poorer  time  resolution  and  viceversa. 

2.1.2  The  WVD 

The  WVD  of  a  signal  s(t)  is  defined  as  the  Fourier  transform  of  the  time-dependent  autocorre¬ 
lation  function 


WVD{t,u) 


exp{—jLot'}dt' 


(2) 


where  s(t  +  t'/2)s*(t  —  t'/ 2)  can  be  seen  as  a  time -dependent  autocorrelation  function.  The 
bilinear  WVD  has  better  joint  time-frequency  resolution  than  any  linear  transform.  It  suffers, 
however,  from  a  problem  of  cross-term  interference,  i.e.,  the  WVD  of  the  sum  of  two  signals 
is  not  the  sum  of  their  individual  WVDs.  If  a  signal  contains  more  than  one  component  in  the 
joint  time-frequency  domain,  its  WVD  will  contain  cross  terms  that  occur  halfway  between 
each  pair  of  auto-terms.  The  magnitude  of  these  oscillatory  cross  terms  can  be  twice  as  large 
as  the  auto-terms.  To  reduce  the  cross-term  interference,  the  filtered  WVD  has  been  used  to 
preserve  the  useful  properties  of  the  time-frequency  transform  with  a  slightly  reduced  time- 
frequency  resolution  and  largely  reduced  cross-term  interference.  The  WVD  with  a  linear 
low-pass  filter  belongs  to  Cohen’s  class. 


2.1.3  Cohen’s  class 


The  general  form  of  Cohen’s  class  is  defined  as 


C(t,u )  =  j  J  s  ^  s*  ^  <j>(t  —  u,  t')  exp{—  ju>t'}dudt'  (3) 
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TF  Transforms 

Formula 

Kernel 

Cohen  Class 

C(f,a;)  =  fff  <1 ><0,  Os  +  0  V  +  0  e  ito-i^^dudTdO 

<W,r) 

Wigner- Ville 

WVD(l,u)  =  J.v  +  0.V  +  0  exp {-juiT}dr 

<fr(0.r)=  1 

Choi-Willianis 

CW(l,ui,ir)  =  jf  Kcw(u -t,T).s  ^ii  +  y'js’  ^ ^  e~*^Tdudr 

M0,t)  =  KcwiO.  t)  =  -  1  exp(-02r2/tr( 

47T-*/2  T*fo 

Pseudo  Wigner 

PWD{t%uj,a)  =  f Ii(t)s  ^  exp{  —  jurr}dT 

<1>(0,t)  =  /i(t)  =  exp{/nr2/2} 

Smooth  Pseudo 
Wigner- Ville  (SPWV) 

SPWD(t,u},a)  =  f  x(l  -  u)PWD(l,ui,n)du 

<J>(0,r)  =  /j(r)  =  exp{yor2/2} 

Cone  Kernel 

CKD(t,u))  =  fff  KCK(t  —  u,t)s  ^  e~^fTdudr 

(  «(r);  |//t|  <  1/2 

4>(0,r)  =  Kck{I,t)  =  < 

1  0;  k/r|>l/2 

Table  1:  Bilinear  Time-Frequency  Transforms 


The  Fourier  transform  of  r),  denoted  as  $(u;,  r),  is  called  the  kernel  function.  It  can  easily 
be  seen  that  if  <h(a;,  r)= 1,  then  o{t,  r)  =  d(t)  and  the  Cohen  class  reduces  to  the  WVD.  Other 
types  of  kernel  functions,  which  lead  to  the  Choi-Williams  distribution,  the  pseudo  Wigner,  the 
smoothed  pseudo  Wigner- Ville  and  the  cone  kernel  distribution  (see  Table  1),  can  be  designed 
to  reduce  the  cross-term  interference  problem  of  the  WVD  [2]. 

2.2  RSO  parameters  of  interest 

The  radar  back-scattering  from  rotating  objects  is  subject  to  Doppler  modulations  that  enable 
us  to  determine  dynamic  properties  of  objects  and  provide  useful  information  about  them.  By 
observing  RSOs  using  a  radar,  useful  parameters  for  their  discrimination  can  be  identified 
as  the  rotation  period  (To),  the  maximum  Doppler  frequency  (.//),„„,  )  and  the  maximum  size 
orthogonal  to  the  LOS(T>jJ.  In  particular  we  can  be  divide  them  into  two  categories: 

•  geometry  dependent  parameters (fDmax,  D±) 

•  geometry  independent  parameters  (To) 

The  former  depend  on  the  acquisition  geometry  whereas  the  latter  does  not. 
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3  Radar  Data  Simulator 

The  objective  of  this  section  is  to  define  the  signal  model  of  the  signal  received  when  the  radar 
observe  a  rotating  object.  In  order  to  do  so  we  exploit  the  geometry  employed  for  Inverse 
Synthetic  Aperture  Radar  (ISAR)  technique.  The  monostatic  ISAR  configuration  will  be  con¬ 
sidered  with  a  brief  review  of  ISAR  theory.  In  particular,  will  be  described  system  geometry 
and  the  received  signal  model. 

3.1  Geometry  and  signal  modelling 

Let  us  consider  the  simple  geometry  shown  in  Figure  3.  The  monostatic  radar  consist  of  a 
colocated  transmitter  and  receiver  (TX/RX)  and  is  located  on  the  origin  of  the  reference  system 
Tc(£i,  £2,  £3).  The  axis  £2  is  aligned  with  the  radar  Line  Of  Sight  (LoS). 


By  considering  the  target  as  a  rigid  body,  the  relative  radar  target  motion  can  be  consid¬ 
ered  as  the  superimposition  of  two  contributions:  a  translation  motion  and  an  angular  motion 
which  are  applied  to  the  center  O  of  the  target.  The  translation  motion  is  denoted  as  the  the 
translational  rotation  vector  S7,r  (t).  Note  that  17, r  (t)  is  always  directed  along  £3 .  The  an¬ 
gular  motions  are  represented  by  the  angular  rotation  vector  ila  (t  ).  The  sum  of  these  two 
rotation  vectors  yields  the  total  angular  rotation  vector  Hr  (f).  The  projection  of  Hr  (f)  on 
the  plane  orthogonal  to  the  LoS  is  called  the  effective  rotation  vector  Hrjj  (t),  it  determines 
the  target  aspect  angle  variation  exploited  in  ISAR  image  reconstruction.  As  demonstrated  in 
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[3],  the  plane  orthogonal  to  Qeff  (f)  is  the  Image  Projection  Plane  (IPP).  The  IPP  is  shown  in 
Figure  4. 


Figure  4:  Image  Projection  Plane 


Then,  a  reference  system  Tx  (x\.  x2, £3)  is  defined,  that  is  embedded  in  the  target  and 
centered  in  O.  The  plane  (xi,x2),  whose  axes  correspond  to  the  cross-range  and  range  coor¬ 
dinates,  coincides  with  the  IPP  The  axis  X2  is  aligned  with  £2  whereas  £3  with  f leff  (t). 

The  reference  system  Tx  is  time  varying  since  the  effective  rotation  vector  changes  its 
orientation  in  time.  In  order  to  define  the  reflectivity  function  of  the  target,  a  new  reference 
system  Ty  (y1,  y2, 2/3)  is  defined  to  be  coincident  with  Tx  at  time  t  —  0.  As  it  will  be  clear  in 
the  following,  the  plane  (2/1 , 2/2)  defines  the  IPP  and  its  axes  correspond  with  the  cross-range 
and  range  coordinates  respectively.  In  particular,  CLeff(t )  can  be  expressed  as: 

ft eff(t )  =  i LoS  X  (ftr(t)  X  i Los)  (4) 

where  iLoS  is  LoS  unit  vector  in  the  Ty  reference  system. 

R0(t)  denotes  the  time  varying  distance  between  the  radar  and  the  reference  point  on  the 
target. 

3.1.1  Received  signal 

For  simplicity,  a  radar  target  is  represented  by  a  set  of  point  scatterers  that  are  primary  reflect¬ 
ing  points  on  the  target.  The  point  scattering  model  simplifies  the  analysis  while  preserving 
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Doppler  features.  In  the  simplified  model,  scatterers  are  assumed  to  be  perfect  reflectors,  re¬ 
flecting  all  the  energy  intercepted.  As  shown  in  Figure  5  the  radar  is  stationary  and  located  at 
the  origin  O  of  the  radar- fixed  coordinate  system  (X,  Y,  Z).  The  target  is  described  in  a  local 
coordinate  system  (x,  y,  z)  that  is  attached  to  the  target  and  has  translation  and  rotation  with 
respect  to  the  radar  coordinates.  The  origin  O  of  the  reference  coordinates  is  assumed  to  be  at 
a  distance  R0  from  the  radar. 


Figure  5:  Geometry  of  the  radar  and  a  target  with  translation  and  rotation 

When  radar  transmits  an  EM  wave  at  a  carrier  frequenxcy  /o,the  received  radar  signal  can 
be  expressed  as  follows: 

SR(t)  =  sT(t )  exp  {— j'47t f0RP(t)}  =  sT(t )  exp  {-j$(t)}  (5) 

where  Rp(t)  —  R0(t)+r(t)  is  the  range  from  the  radar  to  the  scatterer  and  <h(t)  =  — jATrf0Rp(t ) 
is  the  phase  function  at  time  t.  The  modulation  induced  by  rotation  structure  can  be  regarded 
as  a  unique  signature  of  the  target.  Let  us  now  consider  a  simple  case  of  a  rotating  blade  where 
a  =  (3  =  0  and  zq  =  0  as  shown  in  Figure  6.  Where  a  and  (3  are  the  azimuth  and  elevation 
angle  respectively.  Assume  a  scatterer  P  at  (x0,  Do)  on  a  rotor  blade  rotates  about  a  center  point 
O  with  a  rotation  rate  of  Q.  The  distance  from  the  scatterer  to  the  center  point  is  1,  and  the 
distance  between  the  radar  and  the  center  point  is  R0,  which  may  be  a  function  of  time  if  the 
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target  is  moving.  Assume  that  at  t  =  0  the  initial  rotation  angle  of  the  scatterer  in  the  blade  is 
6*o,  and  at  t  the  rotation  angle  becomes  9t  =  90  +  fit  and  the  scatterer  is  rotated  to  (xt ,  yt )• 


Figure  6:  Geometry  of  the  radar  and  rotor  blade  when  a  =  6  =  0 


Because  we  assume  both  the  radar  and  the  rotor  are  on  the  same  plane,  at  time  t  the  range 
from  the  radar  to  the  scatterer  can  be  derived  as 


Rp(t )  =  [R02  +  l2  +  21 R0  sin(6*o  +  Of)]1/2  =  Ro  +  vt  +  l  sin(0o  +  fit)  (6) 

where  v  is  the  radial  velocity  of  point  Q  and  (/ / R0 )2  —y  0  for  far-held.  Thus,  the  radar  received 
signal  becomes 

sR(t)  =  exp  {]  [2 Ttf0t  +  ^ RP(t )]  }  =  exp  {j  [2nf0t  +  $p(t)]}  = 

=  exp  { j ^  [Ro  +  vt]}  ■  exp  {j  [2nf0t  +  sin(6*0  +  fit)]  } 

Let  6*0  =  0,  after  compensating  the  motion  and  removing  the  constant  phase  term  in  Eq.7, 
the  baseband  signal  returned  from  the  scatterer  P  becomes 

l  sin(6*0  +  fit)  \  (8) 

By  integrating  Eq.8  over  the  length  of  the  blade  L,  the  total  baseband  signal  becomes  the 


4l7T 

sB(t)  =  exp  <J  j  — 
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following  : 

sin(Ot  +  0O)  j>  sine 

3.2  RSO  models 

To  simplify  the  problem  we  modelled  targets  as  a  series  of  point  scatterers.  The  point  scatterer 
model  is  used  under  the  assumption  to  observe  objects  having  multiple-structures  with  size 
smaller  than  the  sensor  resolution.  This  model , called  the  point  scatterer  model,  is  widely  used 
in  many  radar  application.  In  fact,  the  electromagnetic  back-scattered  signal  from  a  complex 
target  can  be  thought  of  as  if  it  is  backscattered  from  a  set  of  scattering  centres  on  the  target. 
As  a  result,  the  high  resolution  of  the  sensor  allows  to  such  points  to  be  mapped  as  point- 
scatterers.  The  point- scatterer  model  can  be  related  to  the  electromagnetic  scattering  theory 
through  high  frequency  ray  optics  where  a  set  of  highly  localized  ray  phenomena  can  be  related 
to  a  reflection  or  diffraction  point  on  the  target.  These  points  can  include  specular  reflections 
from  smooth  surfaces,  edge  diffractions  from  edge  and  tips,  as  well  as  multiple  scattering  from 
dihedral  and  trihedral  corner  reflections.  An  8-objects  database  has  been  built  as  shown  in 
Figure  7 


47 T  L  .  „  . 

—  -sm  (fit  +  90) 


(9) 


,  .  .  ,47t  L 

W)  =  Lexp^j-- 
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Figure  7 :  RSO  Database 


4  Feature  Extraction  from  radar  data 

4.1  Introduction 

If  a  target  moves  with  a  constant  rotary  motion  u  =  \ujx,ujy,  ay],  the  instantaneous  frequency 
received  will  be  ft{t)  =  2/0/c[cJ  x  r(t)\h  where  f0  is  the  frequency  of  the  carrier,  r(t)  is 
the  distance  of  one  scatterer  from  the  center  of  rotation  and  h  is  the  radar-target  LOS  (Line 
Of  Sight).  The  amplitude  of  /*(£)  depends  on  the  LOS  direction:  in  fact,  if  n  and  Cj  are 
parallel,/* (t)  would  be  zero  despite  the  target’s  rotation.  The  optimal  case  is  when  the  two 
vectors  are  perpendicular  to  each  other  since  the  /*(£)  is  maximum  and  therefore  its  measure¬ 
ment  is  facilitated.  The  Spectrogram  (S)  of  a  rotating  rigid-body  contains  a  sinusoid  for  each 
scatterer  with: 

•  the  same  Xq,0  =|  io  |  (because  of  the  rigidity  property) 

•  different  amplitudes  fo  depending  on  their  distance  from  the  center  of  rotation  (O) 

•  different  phases  according  to  their  position  with  respect  to  the  radar 
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Thus,  if  we  estimates  one  of  these  sinusoids  we  can  extract  the  information  about  Tq  (and 
so  O)  of  the  observed  RSO.  It  is  clear  that  the  furthest  scatterer  from  O  has  associated  the 
signature  with  the  highest  amplitude  in  S.  By  estimating  such  sinusoid,  we  obtain  also  the 
information  about  the  maximum  Doppler  frequency  ( fomax )•  Moreover,  if  we  select  the  edge 
signatures  ,  no  matter  if  O  is  a  center  of  symmetry,  it  is  possible  to  estimate  the  maximum 
dimension  of  the  Object  Under  Test  (OUT).  NOTE:  with  the  term  maximum  we  do  not  refer 
to  the  real  object  but  its  projection  onto  the  LOS  during  the  Observation  Time(T0&).  The  worst 
scenario  occurs  if  LOS  and  Cj  are  parallel, in  this  case  the  radar  does  not  see  the  rotation  of 
the  RSO  and  it  would  not  be  possible  to  estimate  its  size.  The  Radon  transform  of  a  two- 
dimensional  signal  containing  a  two-dimensional  delta  function  is  a  sinusoidal  pattern  with 
amplitude  corresponding  to  the  distance  of  the  point  from  the  origin  and  the  initial  phase  cor¬ 
responding  to  the  phase  of  the  point  position.Therefore,  it  is  obvious  that  a  sinusoidal  pattern  in 
the  time-frequency  plane  (produced  by  a  time-frequency  representation  of  sinusoidally  mod¬ 
ulated  signal)  will  project  on  a  two-dimensional  delta  in  the  Inverse  Radon  Transform  (IRT) 
[4].  By  using  the  IRT  on  a  multi-component  sinusoidal  signal  the  energy  of  each  sinusoid  is 
concentrated  into  a  point.  The  IRT  of  S  is  evaluated  for  each  element  of  a  vector  of  angles 
a  =  9setSjpt  =  [0°  :  0.1°  :  360°] S'jpl .  Among  all  the  IRT  computed  we  pick  the  one  with  the 
highest  concentration  according  to  the  parameter  M.  This  is  done  because  the  points  in  IRT 
has  the  highest  concentration  when  the  IRT  is  calculated  using  the  same  number  of  periods  Sf 
of  the  Spectrogram  .  The  Concentration  Measure  (M)[6]  calculated  in  corrispondence  of  the 
i-th  Sf  value  is  defined  as  follows: 


mi(p,q)  =  IRT(S,9setSf(i 


(10) 


M;  = 


E„  E0IRi(p^)  +  |iRi(p,?)| 


1  4 


(ID 


EP  q)  +  I1  Ri(p,  q)  I 

The  best  choice  according  to  the  concentration  criterion  is  the  Sf  that  produces  the  mini- 
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mum  value  of  M: 


Thus, 


=  argminM 

Sf 


IR(P)  q)  =  IRT(S,0setS°pt)  =  IRT(S,  a) 


(12) 


(13) 


where  9set  is  a  vector  of  M  equally  spaced  values  between  0  and  2ir  and  IR  is  the  IRT  of 
the  OUT  and  will  contain  a  point  for  every  signature.  Once  we  have  the  IRT  of  S,  we  select  2 
points  ( kj ,  rrij)  j=l,2  that  correspond  to  the  estreme  of  the  OUT  along  the  maximum  observable 
size  direction.  After  the  selection  is  done,  from  the  coordinates  (kj,  rrij)  of  the  selected  point 
(SP)  we  can  estimate: 


•  The  modulation  amplitude  as: 


•  The  modulation  phase  as: 

4>j  =  arctan  (rrij/kj) 


(14) 


(15) 


•  The  modulation  frequency  as: 

,  =a_ 

2t r 

•  The  rotation  period  as: 

A  27rfm  rr\  2?r 

u-^T'Ta-Ji 

with  fs  =  M/Tob 

It  is  then  possible  to  write  the  correspondent  sinusoid: 


(16) 


(17) 


Aj  cos(27t fmt  +  (f>j) 


(18) 


After  a  scaling  transformation  we  obtain  the  correspondent  sinusoid  in  the  Spectrogram 


DISTRIBUTION  A.  Approved  for  public  release:  distribution  unlimited. 


17 


LINEAR  AND  NON-LINEAR  TIME-FREQUENCY  ANALYSIS 
FOR  PARAMETER  ESTIMATION  OF  SPACE  OBJECTS 

domain: 

foL  cos  f y*  +  <f>j\ 

•  the  distance  from  O  to  the  j-th  SP  can  be  written  as 

.  Hi) 

rYj.)  ^  D max 

L  2/o 

where  /0  is  the  carrier  frequency  of  the  radar. 

Finally,  the  estimation  of  the  maximum  dimension  of  the  object  is  given  by 

D±  =  +  Df 

4.2  Features  Extraction  Algorithm 

The  algorithm  has  been  implemented  in  two  versions: 

1 .  Manual  mode 

2.  Automatic  mode 

4.3  Manual  Features  Extraction  Algorithm 

4.3.1  Functional  Block  Scheme 

In  Figure  8  the  flow  chart  of  the  manual  Features  Extraction  Algorithm  mode  is  depicted.  This 
algorithm  mainly  consist  of  5  steps: 

1.  Load  Sr  of  the  RSO  under  test 

2.  Compute  the  Spectrogram  (S)  of  sR 

3.  Compute  the  Inverse  Radon  Transform  of  sR 

4.  Manual  selection  of  the  peaks  in  the  Radon  Domain 
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5.  Selected  Points  (SPs)  -  Features  mapping 


/ - \ 

fomar  max  Doppler 
frequency 


Tn :  rotation  period 
fl :  angular  velocity 


sR(n):  received  signal 

STFT  :  Short-Time  Fourier  Transform 

S(n,k):  Spectrogram 

IRT:  Inverse  2-D  Radon  Transform 

PS:  Points  Selection 

FE:  Features  Extraction 

Figure  8:  Features  Extraction  Alghoritm  Workflow 
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4.3.2  Data  Flow 

In  this  section  the  entire  data  flow  has  been  inspected.  Each  processing  step  will  be  described 
in  term  of  four  items: 

1.  Action 

2.  Motivation 

3.  Input 

4.  Output 


Received  signal  loading 

•  Action:  it  allows  the  user  to  load  the  received  signal  of  a  particular  RSO. 

•  Motivation:  load  the  Received  signal  of  the  n-th  RSO. 

•  Input:  A  number  (n)  that  identifies  a  specific  RSO  contained  in  the  Database. 

•  Output:  Received  radar  signal  (sR) 
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Spectrogram 

•  Action:  The  spectrogram  S  of  the  sR  is  computed  using  the  MATLAB  built-in  function 
Spectrogram.m,  which  corresponds  to  the  squared  modulus  of  the  STFT  of  sR. 

•  Motivation:  S  represents  the  received  signal  sR  in  the  time-frequency  domain  where  its 
Doppler  signature  can  be  observed. 

•  Input:  Received  signal  sR 

•  Output:  Spectrogram  S 


ST(n,k) 


SR(nf— ♦ 

r  i 

STFT 

- » 

/ - 

N2 

^ _ 

Inverse  Radon  Transform 

•  Action:  The  IRT  of  sr  is  computed  using  the  built-in  MATLAB  function  iradon.m  with 
the  optimum  value  of  period  of  the  Spectrogram  Sf  defined  in  Eq.  12. 

•  Motivation:  This  allows  to  concentrate  the  energy  of  each  signature  into  a  point. 

•  Input:  S 

•  Output:  IRT 
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Manual  point  selection 

•  Action:The  user  is  asked  by  the  machine  to  select  2  peaks  during  the  script  execution: 

-  first  the  user  is  asked  to  select  the  furthest  peak  from  the  center  (O)  of  the  image 
IR 

-  next  is  the  user  asked  to  select  the  opposite  peak  to  one  selected  previously  with 
respect  to  ”0”. 

•  Motivation: The  furthest  peak  from  O  correspond  to  the  scatterer  with  the  highest  fo, 
thus  the  highest  distance  from  O 

•  Input:  IRT  of  the  signal 

•  Output:  (ki,mi);(k2,m2) 


SP-features  mapping 

•  Action:  Using  the  relationships  outlined  above  (eq.  14-21)  we  map  the  selected  point  in 
the  radon  domain  into  the  features  space 

•  Motivation:  Estimate  the  features  of  interest 

•  Input:  (k1,m1);(k2,m2) 

•  Output:  fDmax,  D±,  Tn 

4.4  Automatic  Features  Extraction  Algorithm 

This  algorithm  is  a  variation  of  the  version  shown  in  4.3. 

4.4.1  Functional  Block  Scheme 

In  particular,  the  manual  selection  (4.3.2)  is  replaced  by  an  automatic  selection. 

The  automatic  Points  Selection  (PS)  consist  of: 
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1.  A  thresholding  operation 

2.  2D  peak  extraction 

3.  Refinement  of  the  extracted  peaks  using  their  Euclidean  distance  from  O. 

4.  One  of  the  peaks  with  maximum  Euclidean  distance  is  selected,  (k\,m\). 

5.  The  line  ”1”  that  intersects  both  O  and  (fci,mi)is  used  to  fill  a  binary  mask  such  that  the 
line  ”1”  divides  the  mask  in  2  parts,  and  the  part  in  which  (ki,m{)  is  contained  is  set  to 
0.  The  coordinates  of  the  second  peak  ( k2;m2 )  are  selected  after  the  masking  operation. 


Sr(") 


ST(n,k 

S(n,k) 
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Figure  9:  Features  Extraction  Automatic  Algorithm  Workflow 


4.4.2  Automatic  Point  Selection  Flow 

Here,  only  the  PS  processing  step  will  be  described,  since  the  description  of  the  other  parts  of 
the  software  have  already  been  given  in  section  4.3.2. 


Thresholding 
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•  Action:  it  applies  a  threshold  (th)  to  the  IRT  matrix  that  depends  on  the  mean  and 


standard  deviation  of  the  IRT  map. 


th  —  mean2(\IRT\)  +  £  std2(\IRT\),  £  =  3.5  (22) 


where, 


^  N  N 

mean2(\IRT\ )  =  —  EE  |  IRTij 

i=  1  j= 1 


(23) 


std2(|/f?T|) 


\  -  l)5 


N  N 


mean2(|/f?T|)|2 


(24) 


•  Motivation:  it  separates  the  content  related  to  the  object  from  the  contribute  due  to  other 


signal  components  (such  as  noise, artefacts  and  so  on). 


•  Input:  IRT 


•  Output:  thresholded  IRT 


2D  Peaks  extraction 

•  Action:  find  peaks  using  local  maxima 

•  Motivation:  the  coordinates  of  each  peak  in  the  IRT  map  are  used  to  define  their  eu¬ 
clidean  metric 

•  Input:  filtered  IRT 
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•  Output:  a  two  column  matrix  of  the  peaks  coordinates,  peaks=[x;y] 


C  \ 

Peaks 

Extraction 

\ _ / 


Euclidean  distance  calculation 

•  Action:the  Euclidean  distance  of  each  SP  is  computed  with  respect  to  the  center  of 
rotation  O. 

•  Motivation:  the  euclidean  distance  is  the  characteristic  of  the  peak  that  we  need  to 
evaluate  which  peak  is  the  furthest  from  ”0”. 

•  Input:  peaks 

•  Output:  vector  of  euclidean  distances,  eucld  Considering  the  i-th  value  of  the  matrix 
peaks  we  call  Xi  and  yi  its  coordinates.  Given  that,  the  euclidean  distance  of  the  i-th 
peak  is  defined  as: 


eudj\xi,  yi)  =  s J (xi  -  x0)2  +  {yi  -  Vo)2 


(25) 


Peak  selection 

•  Action:  the  peak  with  maximum  euclidean  distance  is  chosen.  If  more  than  one  peak 
has  the  same  distance  from  O,  we  pick  the  one  with  highest  amplitude. 

•  Motivation:  the  peak  with  maximum  euclidean  distance  from  ”0”  corresponds  to  the 
point  with  maximum  Doppler  frequency  fnmax 

•  Input:  eucld,  filtered  IRT 
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•  Output: 


(Ay,  mi)  =  arg  max  (eucld) 
x,y 


(26) 


Peak 

Selection 


Geometrical  masking 

•  ActionTt  is  defined  a  line  1  that  intersects  O  and  (Ay,mi).  Using  line  1  we  fill  a  binary 
mask  such  that  the  line  1  divides  the  mask  in  2  parts  and  the  part  in  which  (fe^mjis 
contained  is  set  to  0.  Then  this  mask  is  used  to  select  a  subset  of  the  peaks  extracted 
previously. 

•  Motivation:  In  order  to  estimate  the  maximum  size  of  the  OUT  we  also  need  to  find 
the  point  of  the  target  that  is  diagonally  opposite  to  the  SP  Our  algorithm  search  for  the 
peak  at  maximum  euclidean  distance  from  O.  Thus, in  case  of  asymmetrical  objects  if 
we  need  to  exclude  the  side  of  the  object  where  the  first  detection  is  located.  Otherwise 
we  would  estimate  the  RSO  size  incorrectly. 

•  Input:  peaks,  filtered  IRT, (Ay,  mi), O 

•  Output:  peaks  subset 
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4.5  Data  Analysis 

After  describing  the  Scenario,  in  section  4.5.1,  the  results  of  the  technique  described  in  section 
4.4  are  be  shown  in  section  4.5.3. 


4.5.1  Scenario  description 

The  parameters  used  to  generate  the  simulated  radar  signal  are  shown  in  table  2.  The  geomet¬ 
rical  and  dynamic  properties  used  to  create  the  RSO  dataset  are  outlined  in  table  3  . 


Simulated  radar  system 


Symbol 

Description 

Value 

fo 

Ca  rrier  frequency 

5  GHz 

B 

Bandwidth 

800  MHz 

T0b 

Observation  Time 

1.4  s 

fs 

Sampl  ing  frequency 

13.2  kHz 

Tr 

Pulse  Repetition  Interval 

75.7  /i  s 

Ti 

Pulse  Width 

1.25  ns 

Table  2:  Simulated  radar  system  parameters 
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RSO  Dataset  Geometrical  and  dynamic  properties 


RSO 

Shape 

Size(x  x  y  x  ^)[m] 

D_i_[m] 

fD, m\Hz] 

Tn[s 

1 

Satellite 

4  x  4  x  22 

22.0907 

6614.1043 

0.35000 

2 

Sphere 

2x2x2 

2 

449.1096 

0.46667 

3 

Cylinder 

2x2x6 

6.3246 

2840.4188 

0.23333 

4 

Cube 

4x4x4 

5.6569 

1270.2739 

0.46667 

5 

Parallelepiped 

2  x  16  x  4 

18.5536 

5369.0052 

0.40000 

6 

Cone 

12  x  6  x  12 

12 

2245.5482 

0.56000 

7 

Disk 

12  x  12  x  0 

12 

4491.0965 

0.28000 

8 

Asymmetric  Satellite 

4  x  4  x  22 

22.114 

4799.8502 

0.70000 

Note:  The  distance  between  the  radar  and  each  RSO  is  Rq=900  Km 


Table  3:  RSO  Dataset  Geometrical  and  dynamic  properties 

4.5.2  Performance  Metrics 

As  performance  metrics  we  consider  the  following: 

1 .  Mean  Percentage  Error  (MPE,%Err) 

The  Mean  Percentage  Error  is  the  mean  along  the  number  of  different  realizations  of 
noise  Nr  of  the  magnitude  of  the  difference  between  the  exact  value  a  of  the  quantity 
being  forecast  and  the  forecast  fr  divided  by  the  magnitude  of  the  exact  value  per  100. 

Nr 


r= 1 


MPE  =  V 


2.  Normalized  Root  Mean  Square  Percentual  Error  (NRMSPE,%NRMSE) 

The  Root  Mean  Square  Error  (RMSE)  is  a  frequently  used  as  measure  of  the  differences 
between  values  (sample  and  population  values)  predicted  by  a  model  or  an  estimator  and 
the  values  actually  observed. 


Normalizing  the  RMSD  facilitates  the  comparison  between  datasets  or  models  with  dif- 
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ferent  scales. 


NRMSE 


RMSE 

a 


(29) 


The  NRMSE  expressed  in  percentual  form  gives  the  Normalized  Root  Mean  Square 
Percentual  Error  (NRMSPE): 


NRMSPE  =  NRMSE  *  100 


(30) 


4.5.3  Estimation  results 

In  table  4  and  table  5  are  listed  the  performance  metrics  (%Err,%NRMSE)  of  the  estimations 
obtained  referring  to  the  scenario  introduced  in  section  4.5.1. 
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RSO  1 

RSO  2 

Tn 

^Dmax 

Di 

Tn 

^Dmax 

Di 

SNR=-5dB 

%Err 

0,412 

0,106 

0,538 

SNR=-5dB 

%Err 

0,298 

6,561 

4,649 

%NRMSE 

0,426 

28,027 

0,558 

%NRMSE 

0,304 

28,027 

5,554 

SNR=0dB 

%Err 

0,154 

0,106 

0,311 

SNR=0dB 

%Err 

0,332 

8,339 

7,687 

%NRMSE 

0,159 

0,154 

0,315 

%NRMSE 

0,345 

8,777 

7,965 

SNR=10dB 

%Err 

0,078 

0,106 

0,243 

SNR=10dB 

%Err 

0,354 

17,246 

14,805 

%NRMSE 

0,079 

0,154 

0,249 

%NRMSE 

0,363 

17,856 

15.432 

SNR=20dB 

%Err 

0,089 

0,106 

0,252 

SNR=20dB 

%Err 

0,361 

18,546 

18,038 

%NRMSE 

0,089 

0,154 

0,257 

%NRMSE 

0,368 

18,835 

18,325 

RSO  3 

RSO  4 

Tn 

^*Dmax 

Di 

Tn 

^max 

Di 

SNR=-5dB 

%Err 

0,064 

15,622 

30,005 

SNR=-5dB 

%Err 

0,021 

3,334 

3,416 

%NRMSE 

0,077 

28,027 

36,871 

%NRMSE 

0,033 

28,027 

3,470 

SNR=0dB 

%Err 

0,069 

0,747 

1.132 

SNR=0dB 

%Err 

0,024 

3,334 

3,424 

%NRMSE 

0,081 

0,880 

1,242 

%NRMSE 

0,037 

3,388 

3,477 

SNR=10dB 

%Err 

0,202 

0,885 

1,252 

SNR=10dB 

%Err 

0,037 

3,334 

3,375 

%NRMSE 

0,209 

1,085 

1,393 

%NRMSE 

0,046 

3,388 

3.427 

SNR=20dB 

%Err 

0,229 

0,679 

0,904 

SNR=20dB 

%Err 

0,046 

3,334 

3,366 

%NRMSE 

0,233 

0,757 

0,971 

%NRMSE 

0,053 

3,388 

3,417 

Table  4:  Performance  results  (continue  in  table  5) 
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RSO  5 

RSO  6 

Tn 

fcmax 

Di 

Tn 

fcmax 

D± 

SNR=-5dB 

%Err 

0,035 

0,313 

0,368 

SNR=-5dB 

%Err 

0,658 

1,749 

1,228 

%NRMSE 

0,044 

28,027 

0,417 

%NRMSE 

0,666 

28,027 

1,239 

SNR=0dB 

%Err 

0,030 

0,313 

0,363 

SNR=0dB 

%Err 

0,247 

1,812 

1,642 

%NRMSE 

0,032 

0,360 

0,410 

%NRMSE 

0,252 

1,839 

1,661 

SNR=10dB 

%Err 

0,031 

0,313 

0,365 

SNR=10dB 

%Err 

0,268 

1,812 

1,621 

%NRMSE 

0,033 

0,360 

0,412 

%NRMSE 

0,269 

1,839 

1,639 

SNR=20dB 

%Err 

0,027 

0,313 

0,360 

SNR=20dB 

%Err 

0,312 

1,812 

1,581 

%NRMSE 

0,027 

0,360 

0,411 

%NRMSE 

0,312 

1,839 

1,596 

RSO  7 

RSO  8 

Tn 

femax 

Di 

Tn 

femax 

D± 

SNR=-5dB 

%Err 

0,055 

10,803 

18,261 

SNR=-5dB 

%Err 

0,551 

0,402 

0,513 

%NRMSE 

0,073 

28,027 

21,196 

%NRMSE 

0,563 

28,027 

0,569 

SNR=0dB 

%Err 

0,114 

0,355 

0,371 

SNR=0dB 

%Err 

0,307 

0,402 

0,509 

%NRMSE 

0,117 

0,423 

0,432 

%NRMSE 

0,318 

0,405 

0,532 

SNR=10dB 

%Err 

0,183 

0,292 

0,405 

SNR=10dB 

%Err 

0,348 

0,402 

0,485 

%NRMSE 

0,184 

0,292 

0,406 

%NRMSE 

0,355 

0,405 

0,518 

SNR=20dB 

%Err 

0,205 

0,292 

0,428 

SNR=20dB 

%Err 

0,368 

0,402 

0,465 

%NRMSE 

0,206 

0,292 

0,428 

%NRMSE 

0,375 

0,405 

0,501 

Table  5:  table  4  continued 


The  worst  behaviour  is  obtained  in  the  case  of  RSO  2.  The  reason  behind  this  is  to  be 
found  in  the  low  value  of  the  Doppler  frequency,  which  is  harder  to  estimate  with  respect  to 
higher  values  due  to  the  Doppler  quantisation  effect.  On  the  other  hand,  the  best  behaviour  is 
produced  by  RSO  1,  whose  Doppler  frequency  is  the  highest. 

By  looking  at  tables  4  and  5,  and  not  considering  the  worst  case  of  RSO  2,  we  note  that  ex¬ 
cept  for  the  case  of  SNR  of  -5dB  the  maximum  %Err  and  %NRMSE  do  not  exceed  (0.4, 3.4, 3.4) 
respectively  for  ,  fDmax  and  Dj_.  Although  this  is  the  case  of  simulated  data,  these  results 
give  some  indications  about  the  accuracy  and  robustness  of  the  parameter  estimator  with  re- 
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spect  to  input  parameters.  A  more  complete  analysis  of  the  algorithm  performance  is  given  in 
terms  on  classification,  as  it  will  be  shown  in  sections  7.2  and  7.3. 
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5  Optical  Features 

In  this  section,  RSO’s  optical  features  of  interest  will  be  defined  and  the  methodology  that  is 
used  to  extract  them  from  measurements  will  be  outlined.  An  introduction  to  photometry  will 
be  presented  in  sec.  5.1,  as  photometry  is  generally  used  to  estimate  some  object  parameters 
such  as  spin  rotation.  Next,  the  software  used  to  simulate  optical  data,  the  reference  frames 
used  and  an  analisys  of  generated  data  will  be  described  in  sec. 5. 3  and  5.4.  A  characterization 
of  the  noise  will  be  discussed  in  sec. 5. 5  and,  finally,  a  summary  of  radar  and  optical  features  is 
provided  in  sec.  5.6. 

5.1  Introduction  to  photometry 

5.1.1  Definitions 

Photometry  is  an  optical  technique,  that  consists  of  measuring  the  flux  of  an  astronomical  ob¬ 
ject’s  electromagnetic  radiation.  At  its  most  basic,  photometry  is  conducted  by  gathering  light 
with  a  telescope  and  then  capturing  and  recording  the  light  energy  with  a  photosensitive  in¬ 
strument  (CCD  is  the  most  used). 

In  the  last  decades,  space  debris  identification  and  classification  has  become  a  critical  task  in 
order  to  prevent  collisions  in  space.  Photometry  has  been  significantly  used  to  extract  geomet¬ 
ric  features  of  this  Resident  Space  Objects  (RSOs)  and  to  classify  them  [14],  [13],  [12]. 

Once  an  optical  observation  is  obtained,  a  preprocessing  operation  is  carried  out,  including  a 
masking  operation  and  a  calibration.  Light  curves  are  extracted  from  the  preprocessed  image, 
which  is  also  called  photometric  signature,  and  a  RSO  classification  can  be  made.  This  curves 
represent  the  magnitude  of  light  flux  as  a  function  of  time  or  as  a  function  of  phase  angle  (PA), 
which  is  the  angle  between  the  light  incident  onto  an  observed  object  and  the  light  reflected 
from  the  object. 

Acording  to  [13]  the  following  classification  can  be  made: 

•  Canonical  Class:  if  there  is  a  zero  PA  peak  with  the  brightness  decreasing  on  either  side 
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of  zero  PA  in  a  symmetric  and  smoothy  monotonic  manner. 

•  A2100  Class:  if  there  is  a  local  minimum  in  brightness  at  zero  PA.  On  either  side,  the 
brightness  increases  smoothly  up  to  about  ±40°  PA.  At  greater  phase  angles,  the  bright¬ 
ness  again  decreases  smoothly.  Most  of  these  signatures  show  general  but  imperfect 
simmetry. 

•  Telstar  Class:  if  there  is  an  underlying  Canonical  signature.  However,  there  are  sec¬ 
ondary  minor  brightness  peaks  (seen  mostly  in  the  blue  portion  of  the  spectrum)  at  about 
±40°  PA. 

•  BSS702C  Class:  if  there  is  an  underlying  Canonical  signature.  However,  there  are  sec¬ 
ondary  major  brightness  peaks  at  about  ±60°  PA.  The  secondary  peaks  are  not  only 
brighter  but  also  broader  than  in  Telstar  Class.  In  this  type,  the  secondary  peaks  are 
known  to  come  from  solar  panel  concentrators. 

Another  classification  is  described  in  [14],  according  to  which  objects  can  be  classified  as: 

•  Rocket  blocks  (R/B) 

•  Abandoned  Satellite 

•  Space  debris 

According  to  [14],  all  Russian  R/B  observed  give  pratically  identical  brightness  curves 
with  the  following  characteristic  indications  for  determining  their  type: 

•  The  shape  of  a  brightness  curve  represents  periodic  structure  containing  two  identical 
peaks  on  main  rotation  period  divided  by  antipeaks. 

•  The  Fourier  spectrum  of  a  brightness  curve  contains  two  main  components  associated  to 
a  rotation  period  T0  and  to  a  second  order  simmetry. 
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•  On  some  Fourier  spectrum  curves,  a  third  component  of  noticeably  smaller  amplitude  is 
visible.  It  appropriates  to  thin  structure  of  peaks  on  period  because  of  both  peaks  and 
antipeaks  have  forked  vertex  in  this  case. 

The  main  characteristic  property  of  nonfunctioning  spacecrafts  (abandoned  satellites)  is 
the  violation  of  their  working  type  of  stabilization  in  orbit.  From  the  point  of  view  of  a  pho¬ 
tometric  signature  it  results  in  additional  periodicity  in  the  variations  of  brightness  for  this 
satellites.  Even  if  the  object  had  earlier  spin  stabilization  the  noticeable  precession  of  its  rota¬ 
tion  axis  occurs  when  objects  change  their  own  state.  This  precession  is  observed  in  variations 
of  satellyte  brightness.  Light  curves  are  then  used  to  determine  the  shape  of  an  object  and  the 
attitude  (functioning  or  abandoned).  Concerning  space  debris,  is  possible  to  select  three  main 
indications,  which  should  be  exhibited  in  photometric  observations  of  this  RSOs: 

•  Rather  small  visible  brightness,  which  is  connected  to  their  small  reflecting  surfaces. 

•  Very  jagged  curve  of  brightness,  connected  to  the  depending  of  the  dissipation  diagram 
of  space  debris  on  a  phase  and  aspect  angle. 

•  Very  poorly  expressed  periodicity  of  brightness  change  that  is  characteristic  for  any  kind 
of  rotated  gauze  structures  and  small  quasi-plane  structures. 


In  [12]  a  mathematical  model  of  light  curves,  based  on  the  object  shape,  is  discussed. 
Particularly,  the  total  brightness  of  a  target  object,  made  up  of  n  facets,  as  seen  by  the  observer 
is  given  by: 


L  =  S(n,  /io,  a)pjdj  (31) 

j= i 

where  S  is  the  scattering  law  (or  Bidirectional  Reflectance  Distribution  Function,  BRDF),  pj 
is  the  albedo  of  the  j-th  facet,  a,j  is  the  area  of  the  j-th  facet,  and  alpha  is  the  solar  phase  angle. 
The  BRDF  model  determine  the  way  in  which  light  is  reflected  off  of  each  plate.  Each  facet 
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has  a  surface  area,  albedo,  and  surface  normal  vector,  which  completely  define  the  shape’s 
reflectance  properties.  The  parameters  //  and  //0  are  equal  to  the  dot  product  of  the  facet 
normal  vector  with  and  ufun,  respectively.  The  solar  phase  angle  (pi)  is  the  observer- 
target-Sun  angle,  defined  by  arccos( •  u^).  The  BRDF  used  in  [12]  is  a  combination  of 
Lommel-Seeliger  (Sj.s)  and  Lambert  ( Sl )  laws,  and  is  constructed  as: 


S(/d,fi0,a)  =  f(a)[SLS(n,  Ho)  +  cSL(fi,Ho)}  = 


f  (o;)p./i  o 


(-L- 

\H  +  pO 


(32) 


where  f(a)  is  the  phase  function: 


f(a)  =  A0exp 


T  ka  T  1 


(33) 


The  parameter  c  is  a  weithing  factor  that  determines  the  contribution  of  the  Lambert  scatter¬ 
ing  law.  The  phase  function  is  determined  by  the  amplitude  and  scale  length  of  the  opposition 
effect  (A0  and  D )  as  well  as  the  slope  of  the  phase  curve  (. k ).  [12]  uses  the  following  param- 
ethers: 


c 

0.1 

Ao 

0.5 

D 

0.1° 

k 

-0.5 

Table  6:  Parameters  used  in  the  BRDF  model  by  [12] 


For  SS  A  purpose  two  other  BRDF  model,  known  as  Ashikhmin-Shirley  and  Cook- Torrance, 
have  been  in  common  use  in  recent  years.  They  have  the  advantage  of  being  fast  empiri¬ 
cal  models  which  incorporate  both  diffuse  and  specular  reflections.  These  models  are  likely 
to  represent  typical  spacecraft  materials  more  accurately  than  the  Lommel-Seeliger/Lambert 
scattering  law. 
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Another  feature  can  be  extracted  fr  om  the  flux  measured,  which  is  the  Albedo-area  prod¬ 
uct.  In  [15]  a  relationship  between  Albedo-area  and  the  flux  is  presented,  albedo- Area  is  the 
product  of  the  albedo  of  the  object  (the  ratio  of  the  electromagnetic  radiation  reflected  by  an 
object  and  the  amount  that  is  incident  on  it)  with  its  area.  While  albedo-Area  (aA)  is  an  in¬ 
trinsic  property  of  an  object,  the  observed  aA  of  it  at  any  given  moment  is  a  projection  of  its 
geometry  of  observation  with  respect  to  the  source  of  illumination  and  the  sensor.  ’’Projected” 
aA  of  an  object,  then  refers  to  its  observed  value  based  on  the  current  geometry  of  observation. 
’’Intrinsic”  aA  of  an  object  refers  to  its  value  independently  of  the  geometry  of  observation. 

Considering  a  surface  AA,  flux  is  defined  as  the  net  sum  of  the  energy  A E  flowing  both 
inward  and  outward  through  AA  over  time  Deltat  as  A  E,  A  A  and  At  tend  to  zero: 


/  =  lim  V  — —  (34) 

AA,AE,At—tO  ^  AAAt 

Flux  is  also  related  to  apparent  magnitude  (unit-less  measure  of  the  brightness  of  a  celestial 
object  relative  to  that  of  another  celestial  object  in  a  certain  electromagnetic  waveband)  as 
follows: 


mobj  -  m Ref  =  -2.5  log10  (35) 

J  REF 

Considering  a  two-facet  model,  in  which  an  object  can  be  decomposed  in  a  panel  compo¬ 
nent  (P)  and  a  body  component  (B),  the  projected  albedo-Area  can  be  writed  as  follows: 


\RS\2\OR\2  mSUN~mRSO  .  ,  .  .  ,  .,  ,  , 

7 r —  — 10  2  5  =  aAp  cos(cc)  cos(0)  +  aAg  cos(r/))  cos('^)  (36) 

I  At/ 1 2 

in  which: 

•  6  is  the  angle  between  the  observer- RSO  vector  and  the  panel  normal  vector,  and  it 
determines  the  panel’s  projected  area  from  an  observer’s  point  of  view. 

•  lo  is  the  angle  between  RS  O-Sun  vector  and  the  panel  normal  vector,  and  it  determines 
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the  panel’s  projected  brightness  per  unit  of  area. 

•  is  the  angle  between  the  RSO-Sun  vector  and  the  RSO-Earth  vector,  and  it  determines 
the  projected  brightness  per  unit  of  area  of  the  body’s  nadir-pointing  face. 

•  r i  is  the  angle  between  the  RSO-Earth  vector  and  the  observer-RSO  vector,  and  it  deter¬ 
mines  the  projected  area  of  the  body’s  nadir-pointing  face  from  the  observer’s  point  of 
view. 

•  \RS\  is  the  RSO-Sun  vector  length  in  meters. 

•  \OR\  is  the  observer-RSO  vector  length  in  meters. 

•  \AU\  is  the  Astronomical  unit  in  meters. 

•  Ap  and  A b  are  the  areas  of  panel  and  body  respectively. 


5.2  Data  generation 


In  order  to  evaluate  the  performances  of  Optical  classifiers  we  want  to  simulate  optical  data 
acquisition.  To  accomplish  this  task  the  Blender  software  will  be  used. 


Figure  10:  Blender  window  with  default  setup 
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Blender  is  a  free  open  source  CGI  software,  which  allow  to  create  3D  models  of  objects 
and  scenes,  render  images  and  animations,  configuring  materials  properties  such  as  surfaces 
reflectivity,  light  sources  and  camera  motions. 

For  our  purpose  we  will  use  Blender  to  generate  the  objects  listed  in  Fig.  7.  In  order  to  obtain 
light  curves,  animation  has  been  generated  in  which  objects  have  both  traslational  and  rota¬ 
tional  motion  and  the  rotation  is  around  a  fixed  axis.  Once  obtained  a  set  of  frames,  each  of 
those  has  been  filtered  with  a  moving-window  filter  and  the  relative  magnitude  of  the  object  is 
taken. 

In  the  following  paragraphs  will  be  discussed:  the  reference  frame  used,  the  generated  data 
analisys  and  will  be  characterized  the  noise  of  optical  measurements. 


Figure  11:  Example  of  optical  measurement  of  satellite  object 


5.3  Reference  frame  and  measurement  units 

In  this  section,  the  reference  frames  used  for  the  simulations  will  be  described.  For  our  pur¬ 
poses,  we  can  define  two  reference  frames:  the  observer’s  reference  frame  and  the  object 
reference  frame.  The  first  one  can  be  used  to  characterize  the  traslational  motion  of  the  object 
and  the  measurements,  the  second  one  can  be  used  to  characterize  the  rotational  motion  of  the 
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object. 

The  observer’s  reference  frame  is  chosen  so  that: 

•  The  origin  is  on  the  observer. 

•  The  x-y  plane  is  identified  by  the  azimutal  plane  where  lies  the  bisector  of  the  Field  of 
View  angle  (FOV). 

•  The  x-axis  coincides  with  the  bisector  of  the  FOV  angle  and  is  positive  pointing  the 
object. 

•  The  z-axis  is  orthogonal  to  the  x-y  plane  and  is  positive  so  that  an  observer  on  the  z-axis 
can  see  the  object  moving  from  right  to  left. 


Figure  12:  Observer  reference  frame 

The  object’s  reference  frame  is  not  used  for  measurement,  but  is  useful  to  determine  the 
rotation  motion  at  start-up  of  the  simulation.  It  can  be  defined  as  follows: 

•  The  origin  is  on  the  center  of  the  object. 

•  The  x-axis  is  on  the  side  of  the  object  and  points  on  its  left. 

•  The  y-axis  is  on  the  front  of  the  object. 
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•  The  z-axis  is  on  the  top  of  the  object. 


Figure  13:  Object  reference  frame 


The  measurement  units  used  in  the  simulations  are  explained  as  follows.  Space  is  mea¬ 
sured,  conventionally,  in  Blender  units  [BU],  which  is  the  unit  of  measurement  used  in  the 
software,  so  velocity  is  measured  in  [BU/s].  Angular  momentum  is  measured  in  [rad/s]. 

5.4  Data  analisys 

As  stated  in  the  last  paragraph,  a  set  of  light  curves  are  generated  using  a  CGI  software.  The 
following  cases  are  analyzed: 

1.  Cube  rotating  around  his  z  axis  at  3.93  [rad/s]. 

2.  Satellite  model  rotating  around  his  z  axis  at  3.93  [rad/s]. 

3.  Satellite  model  rotating  around  his  z  axis  at  3.93  [rad/s]  observed  on  a  different  angle  of 
view  respect  the  precedent  simulation. 

4.  Satellite  model  rotating  around  his  z  axis  at  7.85  [rad/s]. 
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5.  Satellite  model  rotating  around  his  z  axis  at  7.85  [rad/s]  observed  on  a  different  angle  of 
view  respect  the  precedent  simulation. 


x  io'3  Simulation  #1  light  curve 


Figure  14:  Simulation  # 1  light  curve 


Simulated  objects  have  also  a  traslational  motion  that  is  parallel  to  the  observer  y-axis  at 
a  velocity  of  2.5  [BU/s].  Considering  simulation  #1  the  light  curve  is  generated  as  shown 
in  Fig.  14,  which  has  a  spectrum  shown  in  Fig.  15.  To  evaluate  the  spectrum,  a  Chebyshev 
window  has  been  used,  in  order  to  reduce  SLL.  Results  of  simulations  #2  —  #5  are  shown 
in  Fig.  16-23.  A  first  characterization  of  each  object  can  be  made  through  the  light  curve 
spectrm.  Considering  simulation  #2  and  #3  it  appears  that  both  spectra  present  the  same  main 
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peaks.  For  this  reason  observations  generated  through  simulation  #2  and  #3  are  related  to  the 
same  object.  The  same  considerations  can  be  made  on  simulation  #4  and  #5.  A  comparison 
between  spectra  related  to  different  observations  of  the  same  object  are  shown  in  Fig.  24  and 
in  Fig.  25. 


Figure  16:  Simulation  # 2  light  curve 


Simulation  #2  light  curve  spectrum 


Figure  17:  Simulation  # 2  light  curve  spectrum 

These  characteristic  frequencies  are  related  to  the  object’s  rotational  speed  and  to  the  object 
geometry.  In  simulations  #4  and  #5  the  object  has  a  double  speed  respect  to  simulations  #2 
and  #3,  so  the  frequencies  shown  are  double.  In  Fig.  26  a  comparison  between  cube  and 
satellite  model  used  in  simulations  #1  and  #2  is  shown.  Considering  a  threshold  of  -60  [dB] 
on  the  spectrum,  both  have  peaks  at  /  =  ±3.8  [H z],  this  frequency  is  relative  to  the  cube 
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Simulation  #3  light  curve 


Figure  18:  Simulation  A3  light  curve 
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Figure  19:  Simulation  A3  light  curve  spectrum 


Simulation  #4  light  curve 


Figure  20:  Simulation  A4  light  cutye 
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Simulation  #4  light  curve  spectrum 


Figure  2 1 :  Simulation  #4  light  curve  spectrum 


Simulation  #5  light  curve 


Figure  22:  Simulation  #5  light  curve 


Simulation  #5  light  curve  spectrum 


Figure  23:  Simulation  ^5  light  curve  spectrum 
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Light  curve  comparison  of  simulation  #2  and  #3 


Figure  24:  Comparison  of  light  curve  spectrums  of  simulation  #2  and  #3 


Light  curve  comparison  of  simulation  #4  and  #5 


Figure  25:  Comparison  of  light  cur\>e  spectrums  of  simulation  7^ 4  and  A  5 


Light  curve  spectrum  comparison  -  Cube  and  Satellite  model 


Figure  26:  Comparison  of  light  curve  spectrums  of  cube  and  satellite  models 
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rotational  motion  and  is  due  to  the  edge  reflection  of  the  cube.  Satellite  model  have  also  peaks 
at  /  =  ±1.88  [.Hz],  due  to  the  motion  of  the  panels.  Objects  are  so  characterized  by  the  number 
of  peaks  over  a  threshold  and  by  peak’s  frequencies. 

5.5  Noise  characterization 

In  the  last  paragraph,  noiseless  measurements  have  been  considered.  Noise  is  composed  of 
undesirable  signal  components  that  arise  from  various  sources  in  an  electronic  imaging  system. 
The  three  primary  sources  of  noise  in  a  CCD  imaging  system  are: 

•  Photon  noise 

•  Dark  noise 

•  Read  noise 

Photon  noise,  also  known  as  photon-shot  noise,  refers  to  natural  variations  of  the  incident 
photon  flux  on  the  device’s  silicon  layer.  The  number  of  photoelectrons  collected  by  a  CCD 
pixel  exhibits  a  Poisson  distribution  and  has  a  square  root  relationship  between  signal  and 
noise: 


Nph  =  (37) 

where  Nph  is  photon  noise  and  s  is  the  signal.  This  kind  of  noise  cannot  be  reduced  via 
camera  design. 

Dark  noise  arises  from  the  statistical  variation  of  thermally  generated  electrons  within  the 
silicon  conducting  the  CCD.  Dark  current  describes  the  rate  of  generation  of  thermal  electrons 
at  a  given  CCD  temperature.  Dark  noise,  which  exhibits  a  Poisson  distribution,  is  the  square 
root  of  the  number  of  thermal  electrons  generated  within  a  given  exposure  time: 

Ndk  =  y/Dt  exp  (38) 
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where  Ndk  is  the  dark  noise,  D  is  the  dark  current  and  texp  is  the  exposure  time.  High 
performance  CCD  camera  systems  reduce  dark  noise  by  cooling  the  CCD  with  thermoelectric 
coolers  (TECs),  liquid  nitrogen  (LA^)  or  cryogenic  refrigeration  during  camera  operations. 
Read  noise  (Nr)  represents  the  noise  introduced  during  the  process  of  quantifying  the  elec¬ 
tronic  signal  on  the  CCD.  The  major  component  of  the  read  noise  arises  from  the  on-chip 
preamplifier.  Spurious  charge  also  contributes  to  the  overall  read  noise  of  the  imaging  system. 
HCCD  cameras  are  able  to  lower  read  noise  by  emploing  carefully  designed  electronics. 

The  signal-to-noise  ratio  for  a  CCD  camera  can  be  calculated  as  follows: 

SNR  =  I'"2  PQetexp  d\  (39) 

A,  y/(P  +  B)Qetexp  +  Dt  +  N? 

where: 

•  P  is  the  photon  flux  incident  on  the  CCD,  measured  in  [photons/pixel/second]. 

•  B  is  the  background  photo  flux  incident  on  the  CCD,  measured  in  [photon/pixel/second] . 

•  Qe  is  the  quantum  efficiency  of  the  CCD. 

•  texp  is  the  exposure  time,  measured  in  [seconds]. 

•  D  is  the  dark  currenrt,  measured  in  [electrons/pixel/second]. 

•  Nr  is  the  read  noise,  measured  in  [electrons  rms/pixel]. 

•  Ai  and  A2  are  the  extremes  if  the  wavelength  interval  considered. 

5.6  Radar  and  optical  feature  summary 

Concerning  Optical  features,  the  number  of  peaks  in  the  light  curve  spectrum  are  related  to 
the  shape  of  the  observed  object,  such  as  Albedo-area  product.  The  peak’s  frequencies  in 
the  light  curve  spectrum  are  instead  related  to  u,  which,  as  already  mentioned,  is  a  geometry 
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Radar  features 

Optical  features 

Geometry 

dependent 

Geometry 

independent 

Geometry 

dependent 

Geometry 

independent 

«/* Dmax 

D± 

Tn 

Main  peak  Frequency 
Number  of  peaks 
Peak’s  frequencies 
Albedo-area  product 

U) 

Object  attitude 
Object  class 
Object  class 

Table  7:  Radar  and  optical  features 


independent  feature.  Generally  speaking.  Albedo-area  product  is  an  object’s  intrinsic  feature, 
but  its  measurement  depends  on  the  particular  geometry.  As  stated  in  section  5.1,  the  object’s 
attitude  is  not  dependent  on  the  geometry  considered  and  can  be  estimated  directly  by  means 
of  light  curves,  and  so  also  the  object’s  class.  In  order  to  classify  the  object  as  one  of  the 
classes  defined  earlier,  the  main  peak  frequency,  the  number  of  peaks  and  the  secondary  peak’s 
frequencies  will  be  used.  The  parameter  u>  will  not  be  used  because  it  is  a  function  of  both  the 
main  peak  frequency  and  of  the  object’s  shape.  Being  the  latter  an  unknown,  the  estimation 
of  u>  is  not  viable.  The  radar  and  optical  features  are  summarised  in  table  7  for  the  sake  of 
completeness. 
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6  Classifiers  Description 

The  purpose  of  this  section  is  to  describe  the  classifiers  that  will  be  used  with  the  aim  to 
associate  a  detected  object  to  a  certain  debris  class,  or,  more  pecisely,  to  check  if  a  new  de¬ 
tection  belongs  to  the  same  class  of  a  previous  detection  after  one  or  more  orbit  revolutions. 
Among  the  many  types  of  classifiers,  two  have  been  analyzed  in  details,  namely  the  Support 
Vector  Machine  (SVM)  classifier  and  the  Bayesian  classifier.  While  the  two  classifiers  have 
some  differences  in  the  underlying  analytical  aspects,  both  focus  on  building  hypersurfaces 
(or  hyperplanes,  if  the  classes  are  linearly  separable)  to  separate  the  defined  classes  with  the 
minimum  error  possible,  or  minimum  cost  achievable.  To  achieve  that,  the  classifiers  act  on  a 
hyperspace  F  defined  by  the  ensemble  of  the  different  features  that  characterize  each  classes. 
In  the  next  paragraph  a  detailed  description  of  the  classifiers’  principles  are  provided. 

6.1  SVM  Classifier 

Classifying  into  k  classes  is  to  find  a  decision  function  that  assigns  each  element  of  the  dataset, 
i.e.  each  feature  vector  extracted,  to  one  of  the  k  classes.  In  the  case  of  binary  classification 
(k  =  2)  one  of  the  most  effective  approaches  is  represented  by  the  Support  Vector  Machines 
(SVMs)  that  will  be  described  below.  Starting  from  binary  SVMs  is  also  possible  to  obtain 
a  classification  into  k  classes  using  ’’voting  schemes”,  or  extending  their  pristine  formulation. 
The  following  paragraphs  will  describe  these  approaches. 

6.1.1  Binary  SVMs 

SVMs  were  introduced  initially  by  Vapnik  [7]  in  the  80’s  and  represent  a  very  effective  tech¬ 
nique  for  non-linear  classification.  Given  a  set  of  input  vectors  in  an  F-dimcnsional  space 
(where  F  is  the  number  of  features)  the  binary  classification  states  that  each  vector  is  associ¬ 
ated  with  one  of  two  possible  classes.  The  basic  idea  of  the  SVM  is  to  build  a  hyperplane  of 
separation  between  the  two  classes  in  the  F-dimcnsional  feature  space.  The  training  dataset 
consists  of  (xj,  y% }  pairs,  where  the  input  vector  is  x*  e  MF,  and  the  output  is  yt  €  {—1, 1}, 
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where  yt  —  1  identifies  the  class  of  positives  and  yt  —  —1  identifies  the  class  of  negatives. 

6.1.2  linear  SYM,  separable  case 

The  simplest  case  is  that  of  separable  data  separable.  It  is  assumed  that  in  the  F-dimcnsional 
space  exists  a  hyperplane  that  can  separate  all  positives  from  all  negatives,  as  in  Figure  27. 


Separating 


The  hyperplane  is  identified  by  the  points  x  that  satisfy  the  equation  x  •  w  +  b  =  0,  where 
w  is  the  normal  to  the  hyperplane  and  6/||w||  is  the  distance  from  the  origin.  The  margin  is 
defined  as  the  distance  between  the  positive  element  nearest  to  the  hyperplane  of  separation 
and  the  negative  element  nearest  to  the  hyperplane  of  separation.  The  algorithm  searches 
the  hyperplane  of  separation  with  the  largest  margin.  In  mathematical  terms,  the  margin  is 
maximized  when  ||w||2  is  minimized  respecting  the  constraints  for  each  element: 

x.j  •  w  +  b  >  +1  if  y  t  =  +1 
Xj  •  w  +  b  <  —  1  if  yi  —  —  1 

i.e.: 

Vifc-w  +  b)  —  1  >  0.  (41) 


(40a) 

(40b) 
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The  problem  can  be  solved  by  a  Lagrangian  formulation,  which  allows  to  incorporate  the 
constraints  described  above  in  constraints  on  the  Lagrange  multipliers  (a*)  and  also  permits  to 
show  only  the  scalar  products  between  the  training  vectors.  This  allows,  in  the  case  of  non¬ 
linear  problems,  to  use  kernel  functions  that  map  the  data  from  the  original  feature  space  onto 
a  Hilbert  space  where  the  problem  is  linear.  The  Lagrangian  formulation  is  as  follows: 

LP  =  2  llWl|2  •  W  +  6)  +  (42) 

i  i 

then  the  problem  is  to  minimize  Lp  with  respect  to  w  and  b  and  simultaneously  requesting 
that  the  derivatives  with  respect  to  all  the  a:,  (Lagrange  multiplier)  will  be  equal  to  zero,  with 

OLi  >  0. 


6.1.3  linear  SYM,  non  separable  case:  Cortes- Vapnik  formulation 

If  the  data  are  not  separable  the  method  described  above  does  not  allow  to  find  a  solution. 
To  solve  this  problem,  the  constraints  should  be  relaxed  so  that  some  errors  are  tolerated,  or 
some  element  of  the  postitive  class  may  belong  to  the  half-plane  of  the  negative  class  and  vice 
versa.  The  aim  is  always  to  maximize  the  margin  and  therefore  minimize  || w|| 2  but  with  the 
following  constraints: 


x*  •  w  +  b  >  +1  -  &  if  yi  =  +1,  £  >  0  (43a) 

x»  •  w  +  b  <  -1  +  &  if  yi  =  -1,  £  >  0.  (43b) 

The  variables  Q  are  a  measure  of  tolerability.  To  have  a  classification  error,  Q  must  exceed 
unity,  then  JY  Q  represents  an  upper  bound  on  the  training  error.  The  upper  limit  on  the  max¬ 
imum  distance  between  an  element  and  the  correct  separation  hyperplane  determines  the  cost 
associated  with  the  misclassification  of  that  element.  To  consider  the  cost  of  every  classifica¬ 
tion  error,  we  add  a  new  term  to  the  objective  function  to  be  minimized,  ||w||2/2  ,  so  the  new 
function  to  be  minimized  is: 
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i||wl|2  +  c5>  (44) 

i 

with  the  constraints  in  eq.  43,  where  C  is  a  parameter  to  choose:  the  greater  C  is,  and 
the  greater  the  penalty  associated  to  the  errors.  Incorporating  the  minimization  of  the  new 
objective  function  and  constraints,  using  the  Lagrangian  formulation  we  obtain: 

LP  =  ^IIWH2  +  C  X/  &  ~^2ai[yi(xi  ■  W  +  6)  -  1  +  &]  -  T:  l-kti-,  (45) 

i  i  i 

where  C  represents  an  upper  limit  for  the  Lagrange  multipliers  at,  while  ji%  are  the  La¬ 
grange  multipliers  which  force  the  positivity  of  Q.  This  formulation  of  Cortes-Vapnik  (also 
called  cost-oriented)  [8]  makes  sure  that  the  optimization  process  searches  a  compromise  be¬ 
tween  the  margin  and  the  error  of  misclassification  and  it  is  for  this  reason  that  for  this  formu¬ 
lation  will  be  chosen  for  the  problem  under  consideration. 


Figure  28:  Separation  hyperplane  in  features  space  for  non-separable  data 

In  the  event  that  the  decision  function  is  not  linear  in  the  space  of  features,  the  formulation 
described  above  can  be  generalized  using  the  SVM  of  kernel  functions.  Since,  as  noted  above, 
in  the  Lagrangian  only  the  scalar  products  between  the  vectors  of  training  appears,  it  is  possible 
to  repeat  the  same  treatment  also  in  the  nonlinear  case  as  long  as  the  data  are  mapped  into  a 
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Hilbert  space  S  through  a  mapping  :  Rr  — »  S.  In  this  way,  the  training  algorithm  will 
depend  only  from  scalar  products  in  S,  or  by  functions  of  the  form  T(x()  •  $  (  x;y ) .  Then  it 
can  be  chosen  a  kernel  function  such  that  iC(x,,  xj)  =  $ (x, )  •  $(xj)  so  that  in  the  algorithm 
of  training  appear  only  K  and  not  the  mapping  <I>.  In  this  way  the  problem  can  be  solved  in 
a  space  of  bigger  size  in  which  the  problem  is  more  linear,  and  therefore  there  can  be  traced 
back  to  the  case  of  linear  SVM.  Examples  of  kernel  functions  are  the  Gaussian  kernel: 


or  polynomial  kernel 


(46) 


Kfc,  xj)  =  (xj  •  xj  +  l)p  (47) 

where  p  is  the  degree  of  the  polynomial. 

The  formulation  used  for  the  classification  of  targets  is  cost-oriented  SVM  with  Gaussian 
kernel. 

6.1.4  Approaches  to  multi-class  classification 

Starting  from  the  combination  of  multiple  SVM  binary  one  can  obtain  a  classification  into  k 
classes  [9].  There  are  two  possible  approaches: 
one-against-one : 

Are  constructed  k(k  —  l)/2  binary  SVM  classifiers,  each  trained  to  distinguish  between  two 
classes.  Then  each  classifier  is  trained  with  training  data  from  two  classes.  k(k  —  l)/2  hy¬ 
perplanes  are  constructed  which  separate  each  class  from  each  other  class.  The  classification 
on  the  test  set  is  made  by  comparing  the  results  of  all  classifiers  according  to  a  certain  voting 
scheme.  For  example,  the  class  that  receives  the  highest  number  of  predictions  will  be  the  one 
assigned  (max-wins-strategy). 
one-against-all : 

K  binary  SVM  classifiers  are  constructed,  each  trained  to  separate  one  class  from  all  the  others. 
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Each  classifier  is  trained  with  the  training  data  from  all  classes.  Then  we  obtain  k  hyperplanes 
that  separate  each  class  from  the  remaining  k  —  1  classes.  Even  in  this  case  the  classifiers  are 
then  combined  by  comparing  their  decisions  on  a  test  set  according  to  a  voting  scheme.  The 
class  with  the  highest  decision  value  will  be  the  one  assigned. 

There  is  also  a  third  approach  based  on  a  generalized  formulation  of  the  algorithm  SVM: 

k-class  SVM 

This  is  a  generalization  of  SVM  such  that  the  function  of  decision  considers  all  the  k  classes 
together.  The  optimization  problem  then  considers  the  data  from  all  the  classes.  Despite  the 
number  of  classifiers  needed  in  the  case  of  the  approach  one-against-one  is  greater  than  the 
case  of  one-against-all  and  formulation  of  the  k- class  (which  uses  only  one),  each  individual 
optimization  problem  is  solved  in  a  space  of  much  lower  dimensionality,  because  the  training 
data  used  for  each  SVM  only  come  from  two  classes.  Thus  the  convergence  towards  the 
solution  is  easier  in  the  one-against-one  strategy,  and  also  the  calculation  time  is  lesser  [?]. 
For  these  reasons,  the  approach  that  has  been  chosen  for  the  classification  of  targets  is  one- 
against-one. 

6.2  Bayesian  Classifier 

In  addition  to  the  SVM  classifier,  also  the  Bayesian  classifier  has  been  analyzed.  The  Bayesian 
classifier  requires  knowledge  of  some  a-priori  information,  ie  the  probability  of  occurrence  of  a 
given  class  ( c),  i.e.  P(cj),  and  the  probability  of  observing  a  given  vector  of  features  d,  assum¬ 
ing  that  class  has  been  specified,  i.e.  P(d|Q).  If  such  knowledge  was  perfect,  the  Bayesian 
classifier  would  be  optimal,  in  the  sense  that  the  probability  of  error  would  be  minimal.  In 
practice,  however,  we  do  not  know  either  the  frequency  of  the  classes  or  the  conditional  prob¬ 
ability  of  the  features,  and  it  is  therefore  necessary  to  estimate  these  parameters  directly  from 
the  training  data,  thus  leading  to  a  sub-optimal  classifier. 

The  objective  of  the  classification  is  to  determine  the  best  hypothesis  h  within  a  given  space 
of  classes  H,  knowing  the  training  data  d  =  [d \ ,  d2,  d3,  ■  ■  ■  ,  dn\.  In  particular,  it  is  required 
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to  determine  the  most  probable  hypothesis  given  the  training  set  d  plus  an  a-priori  knowledge 
of  the  probabilities  of  the  various  hypotheses  ht  in  the  space  II  (i.e.,  the  probability  of  each 
class).  Obviously  in  the  case  discussed  here,  we  have  no  information  on  the  a-priori  probability 
of  each  class,  and  it  is  therefore  reasonable  to  assign  to  each  of  them  the  same  value. 

Bayes  theorem  states  that 


P(hi  |d) 


P(d|/ii)P(/ii) 

P(d) 


(48) 


In  classification,  our  interest  is  to  know  P(hi\d),  the  a-posteriori  probability,  that  is,  the 
probability  of  having  a  certain  class,  after  looking  at  the  features  d.  Often,  a  common  approach 
is  to  proceed  with  the  maximum-a-posteriori  (MAP)  criterion.  In  particular, 


h*MAP  =  argmaxP(/i|d)  =  argmaxP(d|/i)P(/i)  (49) 

heH  h£H 

where  P(d)  disappears  because  it  does  not  depend  on  h.  In  our  case,  where  we  have 
assumed  all  classes  equally  likely,  we  can  also  cancel  the  term  P(h)  and  we  get  a  maximum 
likelihood  estimation  (Maximum  Likelihood,  ML). 

Under  certain  conditions  it  can  be  shown  [11]  that  every  learning  algorithm  that  minimizes 
the  mean  square  error  between  the  output  hypothesis  and  probabilities  on  the  training  data 
gives  out  a  maximum  likelihood  hypothesis.  The  objective  of  the  Bayesian  classifier  (like  that 
of  any  classifier)  is  to  determine  a  particular  discriminant  function  gi  :  d  — *  M  for  each  ith 
class. 

As  with  all  kind  of  real  observation,  every  feature  d  is  corrupted  by  noise.  An  important 
assumption  is  that  error  has  a  Gaussian  distribution.  The  choice  of  Gaussian  noise,  under  the 
assumption  of  having  many  realization  in  the  training  set,  it  is  reasonable  under  the  central 
limit  theorem.  Otherwise,  the  performance  will  be  slightly  degraded.  A  training  dataset  cor¬ 
rupted  by  Gaussian  noise,  has  the  form  d  =  d  +  n,  where  d  is  the  noise-free  value  and  n 
is  the  noise.  The  simplest  and  most  mathematically  tractable  formulation  is  to  have  lines  as 
objective  functions  g  of  the  lines  (i.e.  the  hyperplanes  in  the  case  of  multidimensional  spaces). 
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Class  assignment 


Figure  29:  Classifier  operation:  discriminant  function  are  extracted  from  features,  comparing  these 
values  a  class  can  be  selected 

The  underlying  assumption  in  this  discussion  is  that  the  covariance  matrices  of  the  various 
classes  are  all  equal.  Otherwise,  the  decision  boundaries  will  be  quadratic  hypersurfaces  with¬ 
out  changing  generality  to  the  approach. 


Defined  the  probability  density  in  the  case  of  multivariate  Gaussian  noise 


Px(x) 


1  r-k(x~u)TZ  1(x~v) 


(50) 


applying  the  Bayes  rule,  the  objective  function  ^(d)  is  no  more  than  the  P(7/,jd),  namely 


9i(d)  =  P{K  |d)  = 


P(d|/ii)P(^) 

P(d) 


=  -Ux-/ii)TSi  Lx-Mi) 


I  ATI 


P(hi) 

P(d) ' 


(51) 


In  this  equation,  d  is  the  vector  of  features,  pt  is  the  average  value  of  features  for  class,  E  is 
the  covariance  matrix  of  the  features  for  each  class,  the  operator  |  •  |  represents  the  determinant 
of  the  matrix,  N  is  the  number  of  classes.  Back  in  the  conditions  of  eq.  49,  the  constant  terms 
P{hf)  can  be  removed,  having  assumed  equiprobable  classes  and  P(d)  does  not  depend  on  //,. 
Taking  the  logarithm  to  obtain  a  formulation  algebraically  equivalent  but  more  manageable  we 
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obtain 

S'; (d) log  =  “(x  -  /ii)TSr1(x  -  Hi)  -  ^  log|(2vr)JVS?;|.  (52) 

This  function  is  the  discriminant  between  classes.  At  this  point  the  areas  of  discrimination 
in  the  multidimensional  space  of  features  F  are  obtained  in  the  following  manner,  as  by  eq. 
49.  A  point  in  space  F  belongs  to  a  class  hi  if  the  objective  function  g,  for  that  class  assumes 
the  largest  value  among  all  the  various  objective  functions. 


Figure  30:  Example  for  a  bidimensional  space;  line  f  represents  the  linear  objective  funtion ,  line  Iiml 
represents  the  linear  function  minimizing  the  mean  squared  error  Points  in  space  represent 
the  noise  corrupted  features. 

In  the  case  of  equal  covariance  matrix  for  each  class,  i.e.  E,;  =  E,  the  problem  can  be 
further  reduced.  The  second  addend  in  eq.  52  can  be  eliminated  because  is  the  same  for  all 
classes,  and  the  discriminant  function  is  reduced  to  the  first  term.  The  expression  in  the  form 
in  eq.  53  is  called  Mahalanobis  distance  (or  its  square).  This  metric  assumes  a  linear  form. 
Performing  the  calculations,  we  have  that 

9i(d)  =  -^E-M  -  2/4E-1d  +  MfE-Vj.  (53) 

Removing  the  first  term,  which  is  a  constant,  we  obtain  a  form  of  the  type  g-fd)  = 
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wfd  +  «’().  In  the  case  of  covariance  matrix  equal  to  the  identity,  this  distance  is  the  (squared) 
Euclidean  distance.  Minimizing  this  term  means  therefore  minimizing  the  mean  square  error. 
Also,  since  this  term  has  a  negative  sign,  one  must  search  for  its  maximum  instead. 

The  Bayesian  classifier  is  natively  designed  work  with  multi-class,  and  is  the  approach  that 
has  been  chosen  in  this  report. 
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7  Performances  Analysis 

The  performances  of  the  classifiers  are  shown  in  this  section.  The  performances  are  evaluated 
and  tested  considering  radar-only  features,  optical-only  features  and  a  fusion  dataset  com¬ 
prising  the  features  obtained  with  both  sensors  types.  The  various  datasets  are  created  with 
different  Signal  to  Noise  ratios  in  order  to  appreciate  variations  of  classification  performances, 
if  any.  The  performances  are  evaluated  separately  for  Bayesian  Classifier  and  SVM  Classifier. 
In  the  case  of  fused  dataset,  the  used  approach  was  the  fusion-of-features  rather  than  fusion- 
of-decisions.  Some  considerations  must  be  made  regarding  the  fusion  datases.  The  SNR  of  the 
two  dataset  are  derived  from  two  different  sensors,  so  it  is  impossible  to  merge  directly  two 
different  types  of  dataset  (optical  vs.  radar).  Nevertheless,  it  is  possible  to  perform  a  compar¬ 
ison  by  using  a  combination  of  the  various  SNRs  in  order  to  obtain  couples  of  optical-radar 
SNRs,  thus  obtaining  a  2D  performance  map  rather  than  a  ID  graph.  In  this  condition  is  it 
possible  to  better  appreciate  the  classification  performances. 

The  classification  approach  is  based  on  a  one-vs-one  strategy.  Given  a  set  of  possible 
classes,  the  objective  of  the  decision  process  is  to  know  if  a  specific  detection  matches  the 
previous  detection  for  a  specific  class,  or  not.  With  this  approach,  we  have  a  binary  decision 
process  (the  present  detection  matches/does  not  match  the  previous  detection.  In  an  operative 
scenario,  if  the  match  is  true,  the  track  is  mantained,  otherwise  a  new  track  is  initialized. 

The  test  was  performed  for  each  of  the  8  classes  versus  the  rest  of  the  possible  observations, 
which  we  consider  as  ’’unknowns”.  The  same  training  and  testing  procedures  were  performed 
with  Bayesian  and  SVM  classifiers.  The  performances  were  measured  in  terms  of  Probabil¬ 
ity  of  Correct  Classification  ( PCc ,  blue  line),  Probability  of  Missed  Classification  ( PM ,  black 
line)  and  Probability  of  False  Alarm  ( Pfa ,  red  line).  Pec  is  calculated  by  counting  the  oc¬ 
currences  for  which  the  desired  target  was  correctly  detected  plus  the  occurrences  for  which 
the  unknows  were  classified  as  unknowns.  PM  is  calculated  by  counting  the  occurrences  for 
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which  the  desired  target  was  classified  as  an  unknown,  whereas  Pfa  is  calculated  by  counting 
the  occurrences  for  which  an  unknown  is  classified  as  the  target.  Based  on  these  considera¬ 
tions,  it  can  happen  that  even  with  a  Pec  of  more  than  80%,  because  we  correctly  classified 
the  unknows,  we  can  obtain  a  100%  Pm  because  we  completely  missed  the  desired  target.  Due 
to  the  number  of  tests  performed,  the  Pec  resolution  is  0.5%,  the  Pm  resolution  is  4%  and  the 
PFa  resolution  is  0.57%. 

7.1  Optical  Only 

The  performances  obtained  with  the  use  of  optical  features  only  are  shown  in  this  section. 

To  create  the  dataset,  the  simulations  were  performed  with  9  different  SNRs.  The  lowest 
SNR  is  20  dB  while  the  highest  is  60  dB.  The  SNR  step  is  5  dB;  For  each  SNR,  25  realizations 
per  class  are  created.  For  the  optical  case,  a  feature  selection  must  be  performed  based  on  their 
characteristics.  More  specifically,  there  are  classes  that  have  a  different  number  of  features 
with  respect  to  the  others,  namely  the  number  of  frequencies  associated  with  peaks  found  in  a 
spectrum.  Depending  on  the  target,  a  number  up  to  12  frequencies  may  be  found  for  a  class, 
whereas  other  classes  may  have  only  one  or  two.  This  differences  in  the  number  of  peaks 
create  issues  with  the  definition  of  the  classifier,  which  expects  a  preset  number  of  features 
at  its  input.  On  the  other  hand,  there  are  some  features  that  do  not  suffer  of  this  problem, 
such  as  the  number  of  peaks  and  the  main  peak  frequency,  which  are  single  valued  features  for 
each  object.  Therefore,  by  using  only  a  reduced  feature  ensemble,  the  inner  separability  of  the 
features  is  not  compromised,  while  facilitating  the  classification  task.  The  features  taken  into 
consideration  are: 

•  dominant  peak  frequency 

•  number  of  peaks  [discrete] 

•  first  peak  frequency. 
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To  perform  the  classifiers’  training,  a  single  SNR  was  chosen,  namely  45  dB,  which  is  an 
intermediate  value  between  the  lowest  and  the  highest.  It  has  to  be  noted  that  the  test  was 
performed  also  with  the  45  dB  set.  The  test  is  not  unnecessary,  because  the  training  algorithm 
does  not  produce  hypercurves  that  perfectly  separates  features,  rather  it  designs  hypercurves 
that  minimize  the  mean  square  error  (MSE).  In  fact,  it  is  possible  to  obtain  some  classification 
error  even  when  the  SNR  of  the  train  set  and  the  SNR  of  the  test  set  match. 

Figures  31  from  38  show  the  results  for  each  of  the  eight  target  types  obtained  with  a 
Bayesian  classifier. 
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Figure  31:  Bayesian,  Target  1 


Figure  32:  Bayesian,  Target  2 
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Figure  33:  Bayesian,  Target  3 


Figure  34:  Bayesian,  Target  4 
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Figure  35 :  Bayesian,  Target  5 


Figure  36:  Bayesian,  Target  6 
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Figure  37:  Bayesian ,  Target  7 


Figure  38:  Bayesian ,  Target  8 
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The  results  that  we  obtain  with  the  Bayesian  classifier  can  be  either  very  good  or  very  bad, 
depending  on  the  target  under  observation.  While  the  performances  are  nearly  independent  of 
the  SNR,  due  to  the  fact  that  the  features  are  quite  separable  in  the  feature  hyperspace,  results 
are  not  very  satisfying.  Only  Target  1,  Target  3,  Target  5  and  Target  7  show  good  perfor¬ 
mances,  save  for  a  not  negligible  Pfa  for  Target  3  and  Target  7.  In  the  other  cases,  the  target 
is  completely  missed.  The  partial  conclusion  is  that  with  optical-only  features  in  combination 
with  a  Bayesian  classifier  the  results  are  not  encouraging  although  some  slight  improvement 
may  be  achieved  by  tuning  the  classifier. 

Figures  39  from  46  show  the  results  for  each  of  the  eight  target  types  obtained  with  the 
SVM  classifier. 


Figure  39:  SVM,  Target  1  Figure  40:  SVM,  Target  2 
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Figure  41:  SVM,  Target  3  Figure  42:  SVM,  Target  4 


Figure  43:  SVM,  Target  5 


Figure  44:  SVM,  Target  6 


Figure  45:  SVM,  Target  7  Figure  46:  SVM,  Target  8 
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The  results  that  we  obtain  using  the  SVM  classifier  with  optical  features  are  better.  The 
various  targets  are  always  correctly  classified  when  the  SNR  is  higher,  while  the  Pm,  while 
being  high  for  the  lowest  SNR,  decreases  with  increasing  SNR.  The  Pfa  is  low  enough,  with 
the  highest  being  8%  for  Target  7  at  the  lowest  SNR. 

In  conclusion,  the  simulations  suggest  that  if  we  have  only  optical  data,  the  SVM  classifier 
produces  better  results  than  the  Bayesian,  indicating  that  the  information  in  the  optical  data 
may  be  sufficient  to  produce  acceptable  results. 

7.2  Radar  Only 

In  this  section  the  performance  obtained  with  the  use  of  radar  features  only  are  listed. 

To  create  the  dataset,  the  simulations  were  performed  with  4  different  SNRs,  which  are  -5 
dB,  0  dB,  10  dB  and  20  dB.  For  each  SNR,  25  realization  per  class  are  performed.  For  the  radar 
case  we  utilized  all  the  three  extracted  features,  therefore  a  feature  selection  was  not  performed. 
This  means  that  the  extracted  features  result  in  a  high  separation  in  the  features  hyperplane, 
and  that  all  the  features  are  present  for  every  class.  The  features  taken  into  consideration  are: 

•  maximum  spectrogram  Doppler  frequency 

•  maximum  object  size,  projected  othogonally  to  the  LOS 

•  rotation  period 

To  perform  the  classifiers’  training,  a  single  SNR  was  chosen,  namely  0  dB.  It  has  to  be 
noted  that  the  testing  was  performed  also  with  the  0  dB  set.  As  in  the  case  of  optical  data  and 
for  the  same  reasons,  the  test  is  not  unnecessary.  So,  also  in  this  case,  it  is  possible  to  obtain 
some  classification  errors  even  when  the  SNR  of  the  train  set  and  the  SNR  of  the  test  set  match. 

Figures  47  from  54  show  the  results  for  each  of  the  eight  target  types  obtained  with 
Bayesian  classifier. 
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Figure  47 :  Bayesian,  Target  1 


Figure  49:  Bayesian,  Target  3 


Figure  5 1 :  Bayesian,  Target  5 


Figure  48:  Bayesian,  Target  2 


Figure  50:  Bayesian,  Target  4 


Figure  52:  Bayesian,  Target  6 
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Figure  53:  Bayesian,  Target  7 


Figure  54:  Bayesian,  Target  8 


The  results  that  we  have  obtained  with  the  Bayesian  classifier  are  comparable  to  the  ones 
obtained  with  the  Bayesian  classifier  using  the  optical  data.  Also  in  this  case,  the  classifier  per¬ 
formances  are  nearly  independent  of  the  SNR  as  the  features  are  quite  separable  in  the  feature 
hyperspace.  Only  Target  1,  Target  2  and  Target  8  show  good  performances.  In  the  other  cases, 
the  target  is  completely  missed.  The  partial  conclusion  in  this  case  too  is  that  we  need  to  try 
other  classifiers  or  feature  combination. 


Figures  55  from  62  show  the  results  for  each  of  the  eight  target  types  obtained  with  the 
SVM  classifier. 
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Figure  55:  SVM ,  Target  1 


Figure  56:  SVM,  Target  2 
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Figure  57:  SVM,  Target  3  Figure  58:  SVM,  Target  4 
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Figure  59:  SVM,  Target  5 


Figure  60:  SVM,  Target  6 


Figure  61:  SVM,  Target  7  Figure  62:  SVM,  Target  8 
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The  results  that  we  have  obtained  using  the  SVM  classifier  with  radar  features  show  an 
improvement  even  when  compared  against  the  results  obtained  in  the  case  of  optical  features. 
The  various  targets  are  always  correctly  classified  when  the  SNR  is  high,  and  we  have  only  a 
pair  of  cases  where  we  have  total  miss  detection  at  the  lowest  SNR. 

In  conclusion,  even  in  this  case  the  simulation  results  suggest  that  if  we  have  only  radar 
data,  when  using  the  SVM  classifier,  the  results  are  encouraging. 

7.3  Fusion 

In  this  section  the  performance  obtained  with  the  fusion  of  optical  and  radar  features  are  shown. 
Even  if  relatively  good  results  can  be  obtained  by  using  the  SVM  classifier,  it  is  interesting  to 
investigate  whether  the  performances  can  improve  further  by  using  a  fused  feature  set,  espe¬ 
cially  when  using  the  Bayesian  classifier. 

To  create  the  dataset,  the  simulations  were  performed  by  coupling  all  optical  and  radar 
SNRs,  therefore  obtaining  a  total  of  36  couples.  This  operation  makes  sense  as  two  different 
sensors  may  not  produce  the  same  SNR  in  their  measurements.  For  each  SNR,  25  realizations 
per  class  are  created.  The  dataset  has  6  features  in  total,  which  are  the  fusion  of  the  selected 
optical  and  radar  features. 

To  perform  the  classifier  training,  we  chose  the  couple  with  SNROPTICal  =  45  dB  and 
SNRradar  =  0  dB.  As  already  stated  this  was  chosen  as  an  intermediate  case  and  the  test 
was  performed  also  with  the  same  SNR  couple. 

The  following  figures  represent  3D  plots  where  the  blue  surface  represents  the  Pec,  the 
black  surface  the  Pm  and  the  red  surface  the  Pfa-  The  two  horizontal  axes  represent  the  Radar 
SNR  and  the  Optical  SNR,  respectively. 

Figures  63  from  70  show  the  results  for  each  of  the  eight  target  types  obtained  with 
Bayesian  classifier. 
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Figure  63:  Bayesian,  Target  1 


Figure  64:  Bayesian,  Target  2 
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Figure  65 :  Bayesian,  Target  3 


Figure  66:  Bayesian,  Target  4 
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Figure  67 :  Bayesian,  Target  5 


Figure  68:  Bayesian,  Target  6 
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Figure  69:  Bayesian,  Target  7 


Figure  70:  Bayesian,  Target  8 


It  is  interesting  to  note  that  the  Bayesian  classifier,  when  using  the  fused  feature  set,  per¬ 
forms  dramatically  better  than  in  the  radar-only  and  optical  only  cases.  The  classification 
performances,  in  fact,  are  nearly  perfect,  save  for  some  small  percentage  points  at  the  lowest 
SNR  couples.  In  this  case,  it  is  clear  that  the  addition  of  three  more  features  to  each  dataset 
creates  a  larger  (with  more  dimensions)  hyperspace  where  it  is  easier  to  separate  the  various 
features,  which  is  then  more  suited  to  be  utilized  with  a  Bayesian  classifier. 


Figures  71  from  78  show  instead  the  results  for  each  of  the  eight  target  types  obtained  with 
the  SVM  classifier. 


Figure  7 1 :  SVM,  Target  1 


Figure  72:  SVM,  Target  2 
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Figure  73:  SVM,  Target  3 


Figure  74:  SVM,  Target  4 


Figure  75:  SVM,  Target  5 


Figure  76:  SVM,  Target  6 


Figure  77:  SVM,  Target  7 


Figure  78:  SVM,  Target  8 
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The  results  obtained  with  the  fusion  dataset  and  SVM  classifier  are  somewhat  worse  with 
respect  to  the  results  obtained  with  the  Bayesian  classifier.  By  inspecting  thoroughly  the  re¬ 
sults,  and  comparing  them  to  the  results  obtained  in  the  two  previous  sections,  it  is  clear  that 
they  replicates  the  SVM  results  for  the  two  kind  of  sensors  separately.  The  surfaces  from  Type 
2,  3,  4,  5,  6  have  a  visible  variation  only  with  growing  optical  SNR  and  only  for  PM.  If  we  take 
a  ID  section  of  these  surfaces,  we  obtain  the  same  results  for  optical  sensor  data  with  SVMs. 
Similarly,  Type  3  and  7  for  radar  data  replicates  the  results  obtained  with  the  ID  analysis.  It 
means  that  the  features  of  fusion  6D  hyperplane  have  the  same  degree  of  separability  of  the 
3D  hyperplanes  generated  by  optical  or  radar  features  only. 

In  conclusion,  while  the  fusion  of  features  using  the  SVM  classifier  does  not  bring  signif¬ 
icant  advantage  (the  performances  are  not  degraded),  the  fusion  brings  a  significant  improve¬ 
ment  of  the  performances  when  using  Bayesian  classifier. 

As  a  final  remark,  based  on  the  simulation  results,  we  can  say  that  if  we  have  observations 
from  a  single  sensor,  the  SVM  classifier  has  better  performances.  If  data  from  both  radar  and 
optical  sensor  are  available,  the  Bayesian  classifier  is  the  optimal  choice.  This  results  must  be 
confirmed  with  real  data  before  such  conclusions  can  be  considered  reliable. 
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8  Conclusions 

The  present  technical  report  has  presented  some  innovative  algorithms  for  estimating  RSO 
parameters.  More  specifically,  by  making  use  of  radar  data  and  time-frequency  analysis,  tech¬ 
niques  have  been  developed  that  are  able  to  extract  and  analyse  the  object’s  Doppler  signature. 
Some  object’s  physical  parameters  have  been  correlated  to  Doppler  signature  parameters  and 
estimation  techniques  have  been  devised  to  estimate  them.  Specifically,  the  object’s  rotation 
speed  and  the  object’s  size  along  a  specific  direction  can  be  estimated  to  provide  additional 
unique  features  to  be  used  to  discriminate  a  specific  object  among  others.  This  concept  has 
been  stressed  as  the  techniques  and  algorithms  developed  during  this  project  aim  to  improve 
data  association  in  tracking  and  initial  orbit  determination  tasks.  The  results  presented  in  this 
technical  report  show  how  accurately  an  object  can  be  discriminated  from  others  when  using 
optical-only,  radar-only  and  joint-optical/radar  data.  The  algorithms  have  been  tested  by  using 
simulation  data.  This  may  represent  a  limitation  in  terms  of  algorithm  validation  as  the  data 
simulated  is  not  close  to  realistic  scenarios.  Nevertheless,  as  a  first  step,  a  controlled  simula¬ 
tion  environment  has  allowed  for  some  initial  considerations  to  be  made.  It  is  important  to  note 
that  the  results  obtained  only  consider  object  characteristics  that  exclude  the  orbital  parame¬ 
ters.  The  features  defined  and  estimated  in  this  work  can  be  merged  with  the  orbital  parameters 
to  define  an  extended  class  of  tracking  algorithms  where  the  data  association  is  aided  by  the 
presence  of  additional  information  about  the  object.  Some  results  have  highlighted  that,  with¬ 
out  using  orbital  parameters,  in  both  the  cases  of  optical-only  and  radar-only  measurements, 
the  Bayesian  classifier  demonstrate  that  such  measurements  are  insufficient  to  discriminate 
objects  in  a  satisfactorily  manner.  Nevertheless,  by  jointly  using  features  obtained  by  process¬ 
ing  optical  and  radar  measurements,  a  high  probability  of  correct  discrimination  is  achieved. 
The  algorithms  developed  need  further  validation  with  real  data  before  a  final  assessment  is 
made.  The  members  who  participated  in  this  project  are  currently  involved  in  another  SST- 
related  project  under  a  govemment-to-government  agreement  between  Italy  and  USA.  Within 
such  project  the  problem  of  data  association  in  RSO  tracking  with  radar  measurements  will  be 
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carried  out  by  exploiting  the  algorithms  proposed  within  this  project. 
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